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FOREWORD 



This thesis is written during an exciting time for particle physics. The Large Hadron 
Collider at CERN (Geneva) has been operational for longer than two years and huge 
amount of data have been taken. At the end of the last year, it was announced that the 
existence or nonexistence of the Higgs boson, the last missing piece of the Standard 
Model of the elementary particles, is going to be confirmed with the collision data of 
2012 by the end of the year. At the end of the year 2012, a two-year shut down period 
will begin to prepare the collider for its full performance goals. 
LHC is the largest and highest energy particle accelerator and collider that has been 
built until now and it is going to be upgraded to achieve even higher luminosities. This 
upgrade project is called the "Super Large Hadron Collider" project and it is scheduled 
between 2013-2018. The studies I present in this thesis are indirectly motivated by 
the sLHC. 

MAMMA (Muon ATLAS MicroMegas Activity) group at CERN conducts R&D 
on micromegas-type detectors for use in the ATLAS detector at the LHC after the 
sLHC luminosity upgrade. In these detectors, high resistivity materials are used as 
anode electrodes for spark protection. Resistivity and dimensions of these electrodes 
are determined through trial and error processes. Main work of this thesis is the 
development of a simulation tool, for understanding the charge spread and discharge 
dynamics on these resistive anode strips. 

This simulation tool is named "Chani". On a rectangular surface, it is possible to 
calculate the amount of charge transport between the different areas of the surface and 
the time needed for the total discharge using Chani. 

Although this thesis is written in the context of gaseous particle detectors, Chani can 
be useful for any research where the understanding of charge transport dynamics is 
important. 
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would like to thank my advisor Prof. Cenap §ahabettin Ozben and co-advisor Prof. 
Serkant Ali ^etin for their support. I would like to thank Assoc. Prof. Veysi Erkcan 
Ozcan for his highly beneficial suggestions and corrections on the thesis and the code. 
I would like to thank Prof. Nazmi Postacioglu for his important contribution about the 
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supported by the Turkish Atomic Energy Agency. 
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SIMULATION STUDIES OF CHARGE TRANSPORT ON RESISTIVE 
STRUCTURES IN GASEOUS IONIZATION DETECTORS 



SUMMARY 

Radiation detection by ionization of the gas atoms as the radiation passes through 
goes back to the end of the 19* century. Since then, gaseous ionization detectors 
and related technologies are being developed extensively. Major discoveries of this 
longer-then-a-century period can be counted as the discoveries of Geiger-Miiller 
counter, proportional counter, multiwire proportional counter and micro pattern 
gaseous detectors. With a Geiger-Miiller counter, it was only possible to count the 
number of particles passing through, however, now with a modern gaseous ionization 
detector, it is possible to measure the energy of the passing particle and reconstruct the 
track followed by the particle with very high precision. 

One of the important technologies of the micropattem era is micromegas (micromesh 
gaseous structure). Micromegas detectors have been shown to have very high position 
and energy resolutions, however, it has also been shown that high spark rates in of 
these detectors makes them impossible to use in high luminosity experiments such 
Large Hadron Collider experiments. 

The Large Hadron Collider is going to be upgraded to have one order higher luminosity 
(from 10-^^ cm"^s"^ to 10-^^ cm"^s"^) within the sLHC (Super Large Hadron Collider) 
project. Due to the increase in the beam luminosity, multiplicity of the particles in 
the collisions is also going to increase by, roughly, one order. For this reason, those 
components of the general-purpose detectors ATLAS and CMS which are not capable 
of handling the expected rate are also going to be upgraded. 

Muon ATLAS Micromegas Activity (MAMMA) group at CERN developes 
micromegas detectors for the upgrade of the ATLAS Muon Spectrometer. In order 
to overcome the spark issues, MAMMA group used high resistivity anode strips in 
their chambers. In these chambers, readout strips were standing below the resistive 
strips, both following the same pattern and isolated electrically, hence, the detector is 
read-out through the capacitively induced charge on the read-out strips. 

Disadvantage of using high resistivity electrodes is that they requires a longer charge 
removal time. As the incoming particle rate, thus the rate of electrons arriving on 
the anode strips, gets higher, there might not be enough time for the charges to 
be removed. Such a case would result a charge-up that reduces the detector gain. 
Therefore, the resistivity and the geometry of the resistive strips must be optimized to 
have a spark-safe detector with high rate capability. With this motivation, a tool for the 
transient simulation of the charge transport on a rectangular surface is developed and 
presented in this thesis. 

This simulation tool is named "Chani" because of the resemblance to the word 
"charge" and is a reference to the name of a character in the science fiction novel 
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Dune. It is coded as a macro to run within the ROOT Data Analysis Framework. 
Principally, Chani divides the surface into the subrectangles and calculates their 
electrical potentials and the currents in between, through predefined time instances. 
As new charges arrive, they disturb the potential distribution of the surface thus cause 
currents. Discharges are also modeled by defining some of the subcells as ground 
connection points with user-defined connector resistances. 

In Chani, electrical potential of each subrectangle is calculated in every time instance 
using the so-called "method of moments". With this technique, an influence matrix is 
first calculated and then used to map the new charge distributions to the new potential 
solutions at each time. Since the matrix calculation is done only once, the technique 

provides fast solutions as desired. 

It is shown in this thesis that Chani is capable of calculating the spread of the 
charge over the surface and the time needed for the total discharge. In conclusion, 
a functional tool for the simulation of transient charge transport and discharge is 
developed. Besides its application in the optimization of the resistive structures in 
micromegas chambers, this tool can find application in any research where the charge 
transport dynamics are crucial. 
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GAZLI iYONlZASYON DEDEKTORLERINDEKI DlRENgLI YAPILARDA 
YUK TA§INIMININ BENZEXiM gALI§MALARI 



OZET 

Radyasyonun gaz atomlarini iyonize etmesine dayali olarak par9aciklarin algilanmasi 
19. yiizyilin sonlanndan bu yana kuUanilan bir tekniktir. O zamandan bugiine, 
gazli iyonizasyon dedektorleri ve ilgili teknolqjiler yaygin olarak geli§tirilmektedir. 
Bu yiiz yildan uzun sure9teki onemli bulu§lar Geiger-MuUer sayaci, orantili sayici 
(proportional counter), goklu-telli orantili sayici (multiwire proportional counter) ve 
mikrodoku gazli dedektorler olarak sayilabilir. Biitiin gazli iyonizasyon dedektorleri, 
gaz dolu bir odadan gegen par§acigin gazi iyonla§tirmasi sonucunda olu§an serbest 
elektronlann elektrodlar araciligiyla toplanmasi ilkesine dayanir ancak farkli gazli 
dedektorlerin, i§lev yeterlilikleri ve gozunurlukleri biiyiik degi§iklikler gostermektedir. 
Omegin, Geiger-MuUer sayaciyla, yalnizca ka§ tane pargacigin dedektorden gegtigini 
bulmak miimkiin olabilirken, modem gazli dedektorlerde gegen parfacigin enerjisini 
ve izledigi yolu yiiksek hassasiyetle belirlemek mumkiindur. 

Mikrodokulu dedektorlerin one gikanlarindan biri 1996'da geli§tirilen micromegas 
(MicroMesh Gaseous Structure, mikroorgii gazli yapi) dedektorleridir. Micromegas 
dedektorlerinin oldukfa yiiksek enerji ve konum foziiniirliigiine sahip oldugu, bu 
dedektorlerin pargacik fizigi deneylerinde ve tibbi goriintiileme uygulamalarinda 
kuUamlabilecegi gosterilmi§tir. Bunlarin yamsira, 2006 yilinda geli§tirilen yigin (bulk) 
micromegas iiretim teknigi ile, dedektoriin iiretim siireci tek agamaya indirilmi§, 
geni§ alanlara sahip micromegas tipi dedektorlerin iiretimi miimkiin kilinmi§tir. Ote 
yandan bu dedektorlerde, kivilcim olu§ma sikliklanmn da yiiksek oldugu goriilmiigtiir. 
Kivilcim, dedektoriin aktif bolgesindeki iyonizasyon yogunlugunun artarak gazm 
iletkenle§mesi, bunun sonucunda da yiiksek gerilim elektronlannin kisa devre olmasini 
takiben ani bir yiik aki§iyla sonu9lanan durumdur. Micromegas dedektoriinde kivilcim 
olugtugunda, orgii elektrodundaki butun yiik anoda akar. Bu da orgii elektrodu 
tekrar yiiklenene kadar (olii zaman) dedektoriin i§levsiz kalmasma sebep olur. Bu 
durum mikromegas dedektoriiniin, Biiyiik Hadron ^arpi§tiricisi deneyleri gibi yiiksek 
pargacik aki§i olan deneylerde kuUammim olanaksiz kilmaktadir. 

Biiyiik Hadron ^arpi§tiricisi'ndaki i§in parlakligi, 2013-2018 yiUan arasmda siirecek 
olan sBH^ (Siiper Biiyiik Hadron ^arpigtiricisi) projesi kapsaminda, §u anki 
hedeflenen diizeyinden bir mertebe (10-^^ cm"^s"^'den 10^^ cm'^s'^'e) yiikseltilecektir. 
I§in parlakligindaki yiikselme, proton - proton garpigmalari sonunda olu§an son 
iiriinlerin sikligmda da kabaca 10 katlik bir arti§a sebep olacaktir. Bu sebeple, BH^ 
iizerindeki genel ama§li dedektorler olan ATLAS ve CMS 'in bu diizeydeki pargacik 
aki§larim yeterli iyilikte algilayamayacak kisimlannda da yiikseltmeler yapilacaktir. 
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CERN biinyesindeki Muon ATLAS Micromegas Activity (MAMMA) grubu, ATLAS 
dedektoriiniin yiikseltilmesinde kuUamlmak uzere Micromegas tipi dedektorler 
geli§tirmektedir. Kivilcim problemlerinin iistesinden gelmek i§in, MAMMA grubu 
yaptigi dedektorlerde anod elektrodlan olarak yiiksek diren9li hatlar kuUanmigtir. Bu 
dedektorlerde yiiksek direngli hatlar ve okuma hatlari elektriksel olarak yalitilmi§ 
bir §ekilde iistiiste durup, ayni deseni takip etmektedir. Bu sayede, kivilcim olu§up 
orgii elektrodu anod elektroduna kisa devre oldugunda dahi, anod elektrodu boyunca 
bulunan yiiksek direng sebebiyle orgii elektrodundaki yiikler bo§almadan kivilcim 
ortadan kalkmakta, ayni zamanda dedektor okuma elektroniginin bagli oldugu okuma 
hatlan da kivilcimla birlikte gelen yiiklerden yalitilmi§ olmaktadir. Direngli anod 
elektrodlarimn kullammi RPC (Resistive Plate Chamber) ve ATLAS dedektoriinde 
kuUamlan TGC (Thin Gap Chamber) gibi farkli gazli dedektor teknolojilerinde 
de goriilmektedir ancak bu dedektorlerde direngli anod bir katman halindedir. 
MAMMA grubunun geligtirdigi dedektorlerde ise, bunlardan farkli olarak yiiksek 
direngli elektrodlar da §eritler halindedir. Bu §ekilde yiik yayilimina bagli olarak 
olu§abilecek gapraz etkile§im (crosstalk) ve bunun sonucunda da olu§abilecek sahte 
sinyalleri engellemek amaflanmigtir. MAMMA grubu, direngli anodlar kuUandiklan 
dedektorlerin yiiksek kazanglarda kivilcima dayamkli olarak galigtigim gostermi§tir. 

Yiiksek direngli anod elektrodu kuUammimn getirdigi sakmca ise, bu teknigin, yiik 

bogalmasi igin gereken zamani uzatmasidir. Dedektore gelen par9acik sikligi, buna 
bagli olarak da anod elektroduna gelen elektronlarin sikligi arttik§a, gelen elektronlarin 
ortamdan uzakla§tirilamadan yeni elektronlarin geldigi bir durum ortaya §ikabilir. 
Boyle bir durumda, dedektoriin gaz kazancmda dii§meye sebep olacak bir "yiiklenme" 
olu§acaktir. Bu sebeple, anod elektrodlarimn boyutlan ve ozdirenci bir eniyileme 
siirecinden gegmelidir. Ancak bu §ekilde, kivilcima duyarsiz ama yiiksek pargacik 
aki§ini algilayabilen dedektorler iiretmek miimkiin olabilir. Bu tezde sunulan ana 
§ali§ma, bu motivasyonla geli§tirilmi§ olan, dikdortgen §eklindeki yiizeylerde yiik aki§i 
ve yiik bo§alma siireglerinin zamana bagli benzetimini yapmaya olanak saglayan bir 
aragtir. 

Geli§tirilen benzetim araci "Chani" olarak isimlendirilmi§tir. Bu isimlendirmenin 
sebebi sozciigiin, Ingilizcede "yiik" anlamma gelen "charge" sozciigiinii gagnftirmasi 
ve de bilim kurgu romani Dune'daki bir karakterin ismi olmasidir. Benzetici 
ROOT veri analizi ortaminda galigmak iizere bir "makro" olarak kodlanmi§tir. Basit 
olarak, Chani dikdortgen §eklindeki yiizeyi alt dikdortgenlere boler; yiizeylerdeki 
yiik dagilimimn sebep oldugu potansiyeli ve alt hiicreler arasmdaki akimlari aynk 
anlar boyunca hesaplar. Hesaplanan akimlarimn iki an arasmdaki zaman araligiyla 
9arpilmasi bir hiicreden diger hiicreye ta§inacak olan toplam yiikii verir. Yiizeye ula§an 
her yeni yiik, yiizey boyuncaki potansiyel dagilimini degi§tirir ve yiizey akimlannm 
olu§masina sebep olur. Yiik bo§alma suregleri ise bazi hiicrelerin "giig kaynagi baglanti 
noktasi" olarak (belirli bir baglanti direnciyle) tammlanmasi sayesinde olur, kaynaga 
giderek ortamdan uzakla§acak olan yiikler, bu noktadaki akimlann benzer §ekilde 
zaman araligiyla garpilmasi sonucunda elde edilir. 

Chani, her alt dikdortgenin potansiyelini her zaman orneginde "momentler 
yontemi" olarak adlandmlan bir teknikle hesaplar. Momentler yontemi, kendine 
e§lenik operatorlerle tanimlanan homojen olmayan diferansiyel denklemlerin, matris 
denklemlerine donii§tiiriilerek, bilgisayar yardimiyla yakla§ik sayisal goziimlerine 
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olanak saglayan tekniklerin genel adidir. Modem bilgisayarlarm yiiksek hesaplama 
giicu sayesinde, bu teknikle yapilan yakla§ik hesaplardan pratik amaflar i9in yeteri 
kadar dogru sonuglar almak mumkiin olmaktadir. Momentler yontemiyle elektrik 
potansiyelin hesaplanmasinda, her bir alt hiicrenin iizerindeki yiik dagilimi sabit 
kabul edilir; bir alt hiicredeki yiikiin diger alt hiicrede olu§turacagi potansiyel ise, 
althiicredeki biitun yiikiin merkezde toplandiginda diger alt hiicrenin merkezinde 
olu§turacagi potansiyel olarak hesaplanir. Alt hiicrenin kendi iizerinde olu§turacagi 
potansiyel ise, bir dikdortgen boyunca sabit bir yiik dagiliminin kendi merkezinde 
olu§turacagi potansiyeldir. Bu varsayimlar 9er9evesinde, benzetimlerin baginda, 
hucrelerdeki yiiklerin diger hiicrelerde sebep olu§turacagi potansiyellerin hesaplan- 
masini saglayan bir tesir matrisi hesaplanir. Yiiksek olgekli hesaplamalar soz konusu 
oldugunda, tesir matrisinin hesaplanmasi uzun siirmektedir. Ancak, tesir matrisi, 
benzetim igin bir defa hesaplandiktan sonra, yiik dagilimlarmm sebep oldugu elektrik 
potansiyelin hesaplanmasi 9arpma ve toplama iglemlerine indirgenmektedir. Bu da 90k 
sayidaki zaman omegi boyunca, her defasinda elektrik potansiyelin hizli bir §ekilde 
hesaplanmasmi miimkiin kilmaktadir. 

Bu tezin giri§ kisminda Avrupa Niikleer Ara§tirmalar Enstitiisii (CERN) kisaca 
tamtilip, Biiyiik Hadron ^arpi§tincisi ve biinyesindeki deneylerden bahsedilmi§tir. 
Daha sonra ATLAS dedektoriiniin kisimlari kisaca a9iklanmi§ ve MAMMA grubunun 
9ali§mlarana deginilmigtir. Giri§ kisminin ardindan gazli dedektorlerin tarihsel 
geli§iminden bahsedilen ikinci kisim gelir. Bu kisimda, mikrodokulu dedektorlere 
gelene kadarki gazli dedektorler, iyonizasyon odasindan ba§lanarak tarihsel olarak 
tamtilmi§ ve genel davram§lan a9iklanmi§tir. Ikinci kismm son boliimiinde, bu 
tezdeki ana 9ali§mamn da motivasyonu olan micromegas tipi dedektorlerin genel 
ozellikleri, kuUanim alanlari ve iiretim teknikleri gibi konulara deginilmi§, son olarak 
da MAMMA grubunun geligtirmekte oldugu diren9li anodlu micromegas dedektorleri 
iizerinde durulmu§tur. Takip eden kisimda, geligtirilen yiik ta§inimi benzetim araci 
a9iklamr. Momentler metodunun matematiksel tamtimimn ardindan, bu teknigin 
elektrik potansiyeli belirleyecek olan Poisson denklemine uygulam§ina yer verilmi§tir. 
Sonraki boliimlerde, benzetim aracinin 9ali§ma prensipleri detayli olarak anlatilmi§ ve 
omek hesaplamalar ve oztutarlilik testleri sunulmu§tur. 

Yer verilen ilk hesaplama omeginde 2 cm' ye 2 cm boyutlannda kare §eklinde, 
ozdirenci 10^ Q./n olan bir yiizeyin merkezi etrafinda kiimelenmi§ 10^ elementer 
yuke e§it miktarda yiik ula§masi durumu ele alinmi§tir. Modellemede yiizey, iki 
dogrultada 15 par9aya boliinmek suretiyle 225 alt yiizeye aynlmi§tir. Benzetim, 200 
ns'ye kar§ilik gelen 2000 zaman adimmda ger9ekle§tirilmi§tir. Merkezdeki yiikiin 
kenarlara ve kogelere yayilma siireci ve yiikiin tamamen bogalmasi i9in gereken sure 
gozlemlenebilmektedir. Yiik bo§alma siirecinin RC devresine benzer §ekilde iistel bir 
karaktere sahip oldugu goriilmii§tiir. Toplam yiikiin yansinin bo§almasi i9in gereken 
siire yakla§ik 20 ns olurken, yiikiin yiizde doksammn bo§almasi i9in gereken siire 
yakla§ik 90 ns olmu§tur. Buna ek olarak, yiizeyin sigasimn ve yiizeyin merkezindeki 
yiiklerin bo§alma siirelerinin benzetimde kuUanilan alt hiicre sayisiyla degi§imi de 
hesaplanarak benzetimlerin yakinsakligi dogrulanmi§tir. 

Sonsuz boyutlu diizlemdeki periyodik yiik dagiliminm zamanla relaksasyonu problemi 
analitik olarak 90ziilebilmektedir. Kare §ekilli diizlemlerdeki periyodik yiik 
dagiliminin zamanla degi§iminin benzetimi Chani araciligiyla yapilmi§ ve elde 
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edilen sonuglar analitik §6zumden beklenen sonugla kar§ila§tirilmi§tir. Benzetim 
parametreleri dogru sefildiginde, sonsuz diizlem iqin beklenen relaksasyon zamani ile 
Chani'den elde edilen zaman arasindaki goreli fark %3.5'ten az olmu§tur. Boylece 
Chani'nin anlamli sonuflar verdigi dogrulanmi§tir. 

Yer verilen U9uncu benzetimde ise MAMMA grubunun geli§tirmi§ oldugu diren9li 
anodlu micromegas dedektor prototiplerindeki direngli §eritlerin boyutlanna ve yiizey 
direncine sahip bir yiizey modellenmi§tir. Yine baglangigtaki bir yiikiin yayilimi ve 
topraga iletimi incelenmi§tir. Yuk bogalma suresi 90k daha fazla oldugu i§in bu 
benzetimde bir oncekinden 1000 kat daha fazla zaman adimi bulunmaktadir. Chani, 
uzun sureli hesaplari hafiza ta§imi olmadan yapabilmektedir. 

Bir diger omekte ise, 18 mm'ye 100 mm boyutlarindaki bir yiizey 14229 altyiizeye 
boliinerek yiik ta§mimi hesabi yapilmi§ ve Chani'nin yiiksek boyutlu matrislerle 
(14229 X 14229) de sorunsuz bir §ekilde galgabildgi gosterilmi§tir. 

Bu §ali§mada gosterilmi§tir ki, geli§tirilmi§ olan benzetim araci Chani ile dikdortgen 
§eklindeki bir yiizeydeki yiiklerin yayilimi ve yiiklerin bogalmasi i§in gereken 
toplam zaman hesaplanabilmektedir. Sonug olarak, zamana bagli yiik taginimi ve 
bo§almasimn modelleyebilen bir benzetim araci geli§tirilmi§tir. Bu araq, MAMMA 
grubunun geli§tirecegi yeni micromegas dedektorlerindeki direngli yapilarin en 
iyilenmesinde kuUanilmanm yanisira, yiik ta§imm dinamiklerinin anla§ilmasimn onem 
arz ettigi her ara§tirmada uygulama bulabilir. 
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1. INTRODUCTION 



Gaseous ionization detectors have been widely used in particle physics experiments 
for over hundred years. There are various types of gas-filled detectors and all of 
them are based on the very same principle: When an energetic enough particle goes 
through a chamber with gas, it ionizes the gas atoms and creates electron - ion 
pairs and these electrons and ions are collected on electrodes to gather the particle 
track information. Although the fundamental principle is the same, different type of 
gaseous detectors differ in many aspects such as geometry, amplification technique or 
production method. 

One major problem that is shared by all high-gain gaseous detectors is spark. High 
gains are desired to achieve high sensitivities, however, when the free electrons in a 
chamber exceed a critical density, high voltage electrodes become effectively shorted 
through the electron cloud. This results in the immediate discharge of the anode 
through the electron cloud, which could be dangerous for the chamber itself, as well 
as the peripheral devices such as read-out electronics or high voltage power supplies. 

High-resistivity materials have been used as anode materials of detectors like RPCs 
[[U or TGCs [2] in order to prevent harms caused by sparks. Recently at CERN, 
MAMMA group implemented the resistive anode idea to micromegas technology and 
built spark-resistant micromegas chambers Instead of a full resistive layer used in 
RPC and TGCs, these micromegas detectors have resistive strips. 

MAMMA group conducts R&D on micromegas detectors for the upcoming luminosity 
upgrade of LHC [4]. In their prototypes, resistive strips with different resistivities and 
geometries are compared to achieve a good performance. The main study in this thesis 
is motivated by the need for a systematic method for the optimization of the resistive 
structures in micromegas chambers. For this purpose, a tool for transient simulation 
of the charge transport on a rectangular surface with finite resistivity is developed and 
applied to a number of cases. 
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In the rest of this chapter, CERN, the Large Hadron Collider and ATLAS experiment 
are briefly introduced; and the MAMMA (Muon ATLAS MicroMegas Activity) 
group's activities are summarized. The second chapter is a historical introduction of 
the gaseous detectors. In third and fourth chapters, the charge transport simulator's 
working principles and applications are presented. Discussions and conclusions are 
given in the following chapters. 

1.1 CERN & The Large Hadron Collider 

The European Organization for Nuclear Research (CERN) in Geneva was established 
in 1954 with 12 founding states. Since then, fundamental particle physics research 
has been conducted at CERN which yielded many important discovery and inventions 
for both science and technology. Significant examples of these are the invention of 
multiwire proportional chamber, the discovery of neutral currents, the discovery of W 
and Z bosons, the invention of the World Wide Web and the observation and capturing 
of antihydrogen atoms |l5l. Currently, there are 20 European member states of CERN 
and numerous non-member states co-operate via special agreements. Turkey holds an 
"observer" status [6], but has recently applied for full membership Q. 

The first particle accelerator at CERN, the 600 MeV Synchrocyclotron (SC) started 
its operation in 1957, and it was followed by The Proton Synchrotron (PS) in 1959, 
the Super Proton Synchrotron (SPS) and the Large Electron-Positron (LEP) collider in 
1989. Besides these major ones, at CERN, there are and have been different particle 
accelerators, which have served various physics research. Since 2008, largest ever-built 
particle accelerator, the Large Hadron Collider collides protons at the highest energies 
achieved to date 

1.1.1 General features and current status of the LHC 

The Large Hadron Collider is a 26.7 km long circular accelerator located under the 
French - Swiss border near Geneva, in a depth varying from 45 m to 175 m according 
to the surface shape. The tunnel that LHC is built in was already there and had 
previously been used for the CERN LEP machine |[8l. 
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Figure 1.1: Schematic layout of LHC Machine and the Experiments 



LHC contains two rings where protons (or ions) travel in opposite directions and eight 
interaction points are available from LEP construction. Four of these eight interaction 
points are used and equipped with the detectors and ground structures: respectively 
at Point 1,2, 5 and 8 ATLAS, ALICE, CMS and LHCb detectors are built [SJ. LHC 
layout with interaction points are illustrated in Figure 1.1. 

Proton beams are injected to the LHC after four stages of acceleration in other CERN 
accelerators: first, they are accelerated to 50 MeV through LINACs, then to 1.4 GeV 
in Proton Synchrotron Booster (PSB) and to 26 GeV in Proton Synchrotron (PS) and 
finally to 450 GeV in Super Proton Synchrotron (SPS) |[8l. This sequence is illustrated 
in Figure 1.2 along with the Pb ion injection chain and former LEP injection chain. 
Having reached 450 GeV, proton beams are than injected to the LHC pipes through 
interaction points 2 and 8 in opposite directions to be accelerated to 7 TeV as shown in 
Figure 1.2. 
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Figure 1.2: LHC Beam Injection Complex IfTOll . 

Main performance goals of the Large Hadron Collider are the peak luminosity of 10^^ 
cm"^s"^ and the center-of-mass collision energy of 14 TeV. This corresponds to the 
energies of 7 TeV per proton beam, necessitating an 8.33 T dipole field which is to be 
achieved by superconducting magnets |[8l. 

After two years of 3.5 TeV per beam runs, recently in April 2012, LHC started its 
4 TeV per proton beam runs producing 8 TeV collisions which is the current world 
record fl ll . LHC is scheduled to begin its 2 years shutdown at the end of 2012 to 
get ready for 6.5 TeV per beam runs in late 2014 and finally to reach its design goal 
of 7 TeV per beam [fTTll . In 2010 runs, peak luminosity of LHC experiments was 
200 X lO^^cm^^s^^ and as per May 2012, peak luminostity of LHC experiments is 
greater than 60 x IQt'^cmT'^s^^ [|12L 

LHC and its experiments provide a powerful tool to probe the last missing piece of 
Standard Model of elementary particles, the Higgs boson, and the many questions 
beyond the Standard Model such as existence of supersymmetry or extra dimensions 
in nature, potential ingredients of dark matter, the cause of CP violation ni3l and 
existence of new fermion families lfT4]| . ALICE (A Large Ion Collider Experiment) is 
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dedicated to heavy ion (Pb-Pb) collisions and LHCb (Large Hadron Collider Beauty) 
is focused on the physics of the bottom quark. CMS (Compact Muon Solenoid) 
and ATLAS (A Toroidal LHC Apparatus) are "general purpose" detectors and many 
questions of today's particle physics including the ones mentioned above are studied 
in these experiments. 

1.1.2 The sLHC Project 

The Super Large Hadron Collider is the upgrade project scheduled for 2013 - 2018 
with the main goal of increasing the LHC peak luminosity from 10^^ cm"^s"^ to 
10^^ cm"^s"^). The upgrade scheme includes both accelerator improvements such 
as replacements in the injector chain and enhancement of the focusing quadrupole 
magnets, and detector development such as the improvement of the ATLAS and CMS 

m. 

1.2 ATLAS Experiment 

ATLAS (A Torroidal LHC Apparatus) at Point 1 (see Figure 1.1) is one of the two 
"general purpose" experiments located at the LHC interaction points. Being a "general 
purpose" detector, ATLAS is expected to be sensitive to all physics produced by the 
collisions at the LHC. 




Figure 1.3: Detector axes, (a) Cartesian coordinates shown from interaction point, (b) 
Polar angle and pseudorapidty (t]). 

Before beginning to describe detector parts, the detector axes should be clearly defined. 
Conventionally, the beam direction is defined as the z-axis. Side-A, the side facing the 
Point 8 (or Geneva), is defined as the positive z and the opposite side of the detector, 
side-C is defined as the negative z. The x-axis points towards the center of the LHC 
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Semiconductor tracker 

Figure 1.4: Computer generated image of the ATLAS detector [[161 . 

ring from the interaction point and the y-axis points upwards. As usual, angle from 
the X-axis is defined as ^ and the angle from z-axis is defined as 6. Finally, the 
pseudorapidity is defined as T] = — lntan(0/2). This coordinate system is illustrated 
in Figure 1.3. 

A computer generated image of the ATLAS detector, where the sub detectors can 
be seen, is given in Figure 1.4. The detector contains three main parts: trackers, 
calorimeters and the muon system. 

1.2.1 Trackers 

When LHC run in its full performance goals, protons will collide once in every 25 ns 
and at each bunch crossing nearly 1000 particles will come out in |t]| < 2.5 region. 
Momenta of these particles are detected with high resolution {opj/Pj = 0.05%) by 
the combination of pixel and microstrip semiconductor trackers (SCT) and drift tubes 
in the Transition Radiation Tracker (TRT) 
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1.2.2 Calorimeters 

Electromagnetic (EM) and hadronic calorimeters of the ATLAS detector aim to 
measure the energies of electrons, photons and jets which emerge from the LHC 
collisions. Additionally, the calorimeter system prevents electromagnetic and hadronic 
showers from reaching to the muon system by absorbing their energies; hence allows 
the muon system to produce clean signals. 

LAr (Liquid Argon) electromagnetic calorimeter covers the region |77| < 3.2 and 
measures electron and photon energies with a resolution of Oe/E = 10%. Barrel and 
end-cap parts of the hadronic calorimeter also cover \t}\ < 3.2 region and measure 
jet energies with Oe/E = 50% uncertainty, and the forward calorimeters cover 3.1 < 
|t]| < 4.9 and make jet energy measurements with Oe/E = 50% uncertainty El. 

1.2.3 The Muon System 

The muon spectrometer of the ATLAS detector is formed by the toroid magnets 
(which give the ATLAS detector its name) and gaseous detectors of four different 
technologies. Momentum resolution of the muon spectrometer is Opj/Pt = 10% Q. 

The magnetic field generated by the toroid magnets bend the muon tracks to help the 
particle identification. Monitored Drift Tubes (MDTs) and Cathode Strip Chambers 
(CSCs) cover the \r}\ < 2.7 and 2.0 < |r]| < 2.7 regions respectively, and perform 
precise measurements of the muon tracks |(2j. In the region where \ r} \ < 1.05 Resistive 
Plate Chambers (RPC's) and in the region 1.05 < |77| < 2.7, Thin Gap Chambers 
provide trigger and second coordinate measurements. 

1.3 Muon ATLAS MicroMegas Activity 

The sLHC luminosity upgrade by the factor of 10 will produce roughly 10 times 
more particle flux on detectors. Currently installed detectors are capable of detecting 
the particle fluxes, at least, up to five times the expected rate from the nominal 
LHC conditions. An upgrade in the ATLAS Muon System is going to be necessary. 
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especially in the forward region (|t]| > 2) in order to maintain a good detector 
performance (H. 

MAMMA (Muon ATLAS MicroMegas Activity) group, established in 2007, proposes 
the large area (approximately Im x 2m) bulk-micromegas technology as a solution with 
the following performance goals: Counting rate capability greater then 20kHz/ crn^, 
single plane detection efficiency greater than or equal to 98%, 100/im spatial resolution 
for the incident angles less than 45°, second coordinate measurement, two-track 
separation at around 1-2 mm distance, approximately 5 ns time resolution, level- 1 
triggering capability and good aging behaviour HI. 

In 2007, MAMMA group built and tested its first medium size (450mm x 350mm) 
bulk-micromegas prototype and this chamber was tested in 2008 with 120 GeV 
pions. In these tests, a spatial resolution of (35 ± 5/im) for 500/im strip pitch and 
(24 ±7)/im for 250/im strip pitch was achieved [fTTll . Later on, in order to overcome 
the spark problems, MAMMA group built micromegas prototypes with resistive 
anodes [j3|. Different from the former examples of resistive anode technique, resistive 
layers in these detectors were segmented like the read-out strips in order to prevent 
charge- spread effects. 
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2. GASEOUS RADIATION DETECTORS 



After observations of photoelectric effect on metals, first experiments in which the 
ionizing effects of ultraviolet light on air is performed by Lenard in late 19th century. 
Lenard also measured that "negative ions" of gas moved 2000 times faster than the 
positive ions [1181 . Besides being an early clue for the existence of electrons, this 
velocity asymmetry between negative and positive ions is an essential factor for the 
pulse shape on the output of gaseous detectors. 

Schematic drawing for a simplistic gaseous radiation detector is given in Figure 2.1. 
When any type of radiation that is energetic enough to ionize the gas confined between 
the electrodes of different polarity passes through, emergent electrons and ions are 
drifted towards the anode and cathode respectively hence produce an electrical current. 



When the primary electrons that are freed by the radiation are energetic enough to 
ionize the gas molecules again, secondary electron - ion pairs are produced. The 
bias voltage applied on electrodes can accelerate the electrons such that they can 
collide with the gas atoms and make new electron - ion pairs. As this process is 
repeated, an output signal proportional to the energy of the passing particle is produced 
via the so-called "avalanche multiplication". Further increase in the bias voltage 




^ Ionizing 
Radiation 



Figure 2.1: Schematic drawing of an ionization chamber. 
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Figure 2.2: Dependence of the output signal shape on the bias voltage and energy of 
the passing particles for particles with energies Ei and E2 (£'1 > £'2)- 

may result in a situation where every ionization produces sufficient number of free 
electrons that could make the gas effectively a conductor resulting total discharge of 
anode. This region of operation is called "Geiger - Muller region" and is suitable 
for "counting" the number of particles that pass through the gas chamber without 
knowing their energies. The operating region where no secondary ionization occurs 
is called "ion chamber region" and the region where the bias voltage is so low that 
some of the primary electron - ion pairs recombine and do not reach to the electrodes 
is called "recombination region". When the bias voltage is enough for the avalanche 
multiplication to occur, but less than the Geiger-Miiller limit, output signal strength 
becomes proportional to the passing particle energy; and this operating region is called 
"proportional region". Dependence of the output signal strength on the passing particle 
energy and the bias voltage in gaseous detectors for different regions of operation is 
illustrated in Figure 2.2 |fT9l . 

In principle, any type of gas can be used in gas filled detectors as radiation can ionize 
any type of gas. In fact, as mentioned above, ionizing properties of radiation on gases 
were first observed on air. However, there is a certain threshold potential, depending 
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on the type and pressure of the gas, below which the avalanche multiplication cannot 
occur. Since this threshold potential is lower for them, noble gases are commonly used 
as active gases in gaseous detectors. Most common example is argon because of its low 
cost. Another thing that should be taken into consideration is that, an ionized (excited) 
gas atom eventually recombines with an electron and emits an ultraviolet photon in this 
process. This emitted photon can hit the walls of the cathode electrode of the detector 
and free an electron via photoelectric effect and this new free electron can also get 
accelerated with the bias voltage and produce a "fake" signal. In order to overcome 
this issue, a "quenching" gas with a high absorbance for ultraviolet photons is used 
in the gas mixture of detector. Polyatomic gases such as C//4 or CO2 are common 
examples of quenching gases [fT9l . 

2.1 Early Examples 

Primitive examples of gas filled radiation detectors are ionization chambers, 
proportional counters and Geiger - Miiller counters which work in the three different 
region of operation that are explained above (see Figure 2.2). Essential differences 
between these detectors are their region of operation and output signal characteristics. 

Figure 2.1 can be thought of an example of an ionization chamber with parallel plate 
geometry. Since only the electrons and ions from first ionization contribute to the 
current; output signal in this kind of detectors is proportional to the "intensity" of the 
incoming radiation. 

An electrical counter for the alpha particles were first introduced by Rutherford and 
Geiger in 1908 [20]. Later on, in 1928, Geiger and Miiller developed this device that 
is now called a "Geiger - Miiller counter" BTTl l. A schematic view of a Geiger - 
Miiller counter is shown in Figure 2.3. Since any ionization of the gas in a Geiger 
- Miiller tube results in a very dense avalanche due to very high bias current, a high 
current flows through anode to cathode by the help of ionization electrons. Since a 
sharp pulse is created for each passing particle, it is possible to count the number of 
passing particles with a Geiger counter; however, there is no information about the 
energies of the incoming particles (see the "Geiger - Miiller region" on Figure 2.2). 
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Figure 2.3: Geiger-Miiller Counter. 



The proportional counter was invented in the late 1940s [[22| . This device, in principle, 
looks like the Geiger - Miiller counter, however it is operated with lower a bias voltage 
such that the output signal amplitude is proportional to the number of electrons from 
the initial ionization cluster. Since the number of electrons from the initial ionization is 
proportional to the energy of the travelling particle, it is possible to measure the energy 
of passing radiation with the proportional counter. 

2.2 Multiwire Proportional Chamber 

A revolutionary development in the area of gaseous radiation detectors was the 
invention of the Multiwire Proportional Chamber (MWPC) by George Charpak in 
1968 ll23l . For this invention, Georges Charpak was awarded with the Nobel Prize 
in Physics in 1992. 

Multiwire proportional chamber, sketched on Figure 2.4, can be thought of an array of 
proportional counters operating in parallel. It contains a number of anode wires that are 
separated by a distance on the order of 1 mm and typically several centimeters above 
and below them, cathode planes (in the original production of Charpak, meshes are 
used as cathode electrodes) are present. Besides providing energy information of the 
passing particle, important feature of the MWPC is that it can also give the information 
about the track of the particle on the axis perpendicular to the anode wires. Several 
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Figure 2.4: Multiwire Proportional Chamber. In reality, all of the anode wires are 
connected to amplifiers, only one is shown in this figure for the sake of 
clarity. 

instances of ionization could occur during the passage of one particle which can result 
in avalanches on more than one anode wire (as shown in figure) and the output signal 
time and shape depends on the time and position (distance from the particular anode 
wire) of the initial ionization. These pieces of information could be combined to obtain 
the time and the position of the initial ionizations hence the track of the particle. 

MWPCs are commonly used in particle physics experiments. One example is Thin 
Gap Chambers used for the triggering and the second coordinate measurement in the 
ATLAS muon system. In addition to the elements shown in Figure 2.4, TGCs have 
high resistivity (1 MQ./n and 0.5 MQ./D) layers on their cathode planes for spark 
protection 

Although it is extensively developed and used, MWPCs suffered from several 
drawbacks; these were, mainly, difficulties in positioning the anode wires closer than a 
few milimeters, mechanical instabilities due to the electrostatic repulsion between the 
anode wires and slow removal of the positive ions from the active region of the detector 
hence the disturbance of the electric field due to the diffusing ions [20]. In order 
to overcome these limitations, new types of gaseous detectors called "Micropattern 
gaseous detectors" have been developed since the late 1980s. 
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2.3 Micropattern Gaseous Detectors 

The era of the Micropattern Gaseous Detectors (MPGDs) began with the invention 
of the Microstrip Gas Chamber (MSGC) by Oed in 1988 El. This device had thin 
(10 jUm) anode and thick (90 jUm) cathode strips adjacently patterned on an insulating 
substrate and a back-plane electrode below the substrate. The structure is illustrated 
in Figure 2.5. 200 /im strip pitch, one order of magnitude lower than the MWPCs is 
accurately achieved via photolithography production technique. 




10 nm 

back plane 

Figure 2.5: Schematic drawing of the Microstrip Gas Chamber. Field lines around the 
anode are roughly sketched. 



Thin anode electrodes produce an intense electric field such that the avalanche 
multiplication occurs near the anode and thick cathode electrodes on both sides of the 
anode makes the removal of residual ions easy in MSGCs. For a proper operation, 
voltage on the back plane is arranged such that the field lines do not end on the 
insulating substrate. However, even though such an arrangement is done, some of 
the ions end up on the dielectric substrate and cannot be removed when the detector 
is operated at high gains f25\ . This results in field-disturbing charge buildup and 
gain-reducing aging effects which are major issues with MSGCs. With their high 
read-out granularity and fast operation, MSGCs are used in several experiments such 
as HERA-B at DESY and CMS at CERN. 

There are a few alternate micropattern structures, similar to MSGC, that will not be 
explained here. These are microgap chambers Il26ll27l . smallgap chambers [|28ll29ll 
and microdot chambers llBOllBTI . A detailed review on the micropattern detectors can 
be found in the reference ll32l . 
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Since the recent developments in micromegas detectors is the motivation of the work 
in this thesis, the detailed explanation of micromegas and related technologies is given 
in a seperate section. 

One last technology that should be mentioned under the micropattem gaseous detectors 
is the so-called Gas Electron Multiplier (GEM). Introduced by Sauli in 1996 [|33]| . 
GEM is a three layer metal (copper, 18 /im thick) - insulator (polymer, 25 /im thick) 
- metal (copper, 18 jUm thick) structure with through-holes separated by a 100 jlm 
pitch. Applying a 200V potential difference between metal layers, an electric field 
of 40 kV/cm is achieved in the center of GEM holes. Using this technique in the 
parallel-plate geometries it is possible to confine the avalanche multiplication to the 
small region in the GEM holes hence prevent the ion diffusion as much as possible. 

2.4 Micromesh Gaseous Structure 

Micromegas (Micromesh Gaseous Structure) was introduced in 1996 as a new detector 
concept with which the problems with the MSGCs were essentially resolved Il34l . The 
structure (illustrated in Figure 2.6) had 150 jUm anode strips separated by a 200 jUm 
pitch, deposited on a Kapton substrate via metal deposition techniques. 100 /im above 
the anode strips, "the micromesh", a grid of 3 /im thick metals with 17 x 17 /im 
openings stands on the quartz spacers. 3 mm above the micromesh a bigger mesh is 
placed as the drift electrode on top. 

In the micromegas detector, anode potential was set to volts. Voltages of the mesh 
and drift electrodes were set to the negative values such that the electric field between 
mesh and the anode strips was 100 kV/cm and between drift and mesh were 1 kV/cm. 
As a result of this arrangement, avalanche multiplication of the electrons was only 
happening in the small amplification gap between the mesh and the anode strips. The 
structure provided high performance parameters: Signals faster than 1 ns and gains up 
to 105 are reported in [l34ll . 

The challenging process in the construction of a micromegas structure was gluing 
the mesh on the supporting pillars which were placed on the anode plane via 
photolithography. Flatness of mesh, which in this process is determined by the 
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Figure 2.6: Schematic view of the micromegas detector. 

accuracy of the pillars, is essential for the homogeneity of the gain along the detector 
area. In 2006, a new production method, adopted from the Printed Circuit Board (PCB) 
technology, was introduced with the name "bulk micromegas" ||35]| . In this technique, 
anode plane with the copper readout strips is covered by a photoresistive material and 
then the woven wire mesh (used instead of the electroformed micromesh) is placed 
on top. Then the photoresistive materials were etched via photolithography to form 
the pillars. This technique reduced the number of production steps to one and made it 
possible to design large-area micromegas detectors for larger scale applications. 

Several high energy physics experiments are equipped with the micromegas detectors 
such as COMPASS [|36ti4T1l and CAST ||42ti50l experiments at CERN and T2K [|511- 
[531 experiment in Japan. Besides its use in particle physics, Micromegas detectors are 
also proposed to be used in medical and industrial imaging II54[|55II . 



It is reported for the micropattern gaseous detectors that when the number of electrons 
in an avalanche goes beyond a value between 10^ and 10^, the so-called "Raether 
limit" [|56ll . secondary avalanches exceed the electron cloud both from the front and 
the back; resulting in sparks [[571. For a heavily ionizing radiation which would create 
approximately 10^ electron-ion pairs per cm, this limit is reached with a detector gain 
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Figure 2.7: Micromegas structure with the resistive strips. View from the cross section 
along the axis perpendicular to the readout strips. Drift electrode which is 
not shown here stands 4 or 5 mm above the mesh. 

of 10^ - 10^. It is shown that in the case of a spark, all the charge on the mesh of a 
micromegas detector is coursed down on the anode strips [| 58l . A detailed study on 
the sparking in micromegas detectors showed that the spark rate in heavily ionizing 
environment was too big to get micromegas detectors to operate at high particle rates 
such as those encountered under LHC conditions without making any improvements 
on the basic micromegas structure [|59ll . 

In order to build a micromegas detector that would meet the conditions of the Super 
Large Hadron Collider, MAMMA group modified the micromegas structure by placing 
high-resistivity strips above the read-out strips and an insulator layer in between [Q. 
This structure is depicted in Figure 2.7. Both resistive and readout strips were 150 jUm 
wide and 100 mm long with a strip pitch of 250 jlm. 

MAMMA group presented 3 prototypes of their resistive anode bulk-micromegas 
detectors which were same in geometry, but had resistive strips with different 
resistivities. Both resistive and readout strips were at potential and mesh and the 
drift electrodes were connected to the negative high voltages as usual. They compared 
their results with another prototype of same dimensions but without the resistive strips. 

With the high flux experiments, MAMMA group showed that their micromegas 
chambers with resistive strips were stable against the sparks for gas (avalanche) gains 
up to 20,000. High voltage on the mesh was stable and spark currents were lesser than 
a few 100 nA which was more than 1 order of magnitude lower than a chamber without 
resistive strips under the same conditions. 
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MAMMA group recently started to develop large area micromegas chambers with 
longer resistive strips. As the length of the resistive strip increases, its capacitance and 
maximum resistance also increase, thus these parameters should carefully be arranged 
such that possible "charge-up"s on resistive strips due to high rates could be avoided. In 
the next chapter, a simulation tool developed for this purpose is going to be explained. 
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3. SIMULATION OF CHARGE TRANSPORT AND DISCHARGE 



Main study of this thesis is the transient simulation of charge transport on a rectangular 
surface. As explained in the previous chapter, motivation of this study is to optimize 
the resistive structures in the micromegas chambers that is being developed by the 
MAMMA group at CERN. 

Simulation of charge transport is carried out by the solution of the Poisson equation 
through the surface of interest at each time step and calculation of the currents 
between small subcells of the surface. Since an analytical solution of the Poisson 
equation does not exist for a rectangular surface of finite dimensions, an approximate 
method, so called "moments method" is used for the field calculations. In the next 
subsection, general mathematical description of method of moments and, in particular, 
its application to the Poisson equation is introduced. In the following section, general 
working principles of the simulator is explained and some calculations are presented. 

3.1 Electrostatics via Moment Methods 

As it is going to be explained, method of moments is the collection of methods in 
which the self-adjoint linear operators are written in the matrix forms. For the most 
cases, these matrices are infinite dimensional; however, large, but finite, dimensional 
estimations provides accurate enough solutions for practical purposes. General theory 
of linear function spaces are not going to be explained here, but can be found in many 
text books on applied mathematics for physics and engineering problems II60[|61II . 

Moment methods have a wide range of application in electrostatics, electrodynamics, 
microwaves, antennas and network problems. One can find the detailed text covering 
application of moment methods to all these subjects in referred monograph [1621 . Here, 
the general formulation of the problem and its application to the Poisson equation is 
presented as it is in the reference [l62ll . 
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3.1.1 Mathematical description 

Consider the following inhomogeneous equation of the self-adjoint operator L: 

Lf = g (3.1) 

In this equation, the function g is called the source and function / is called the field 
or response. Once the matrix form of the L is obtained and one of the functions / or 
g is known, it becomes possible to calculate the unknown function directly or after a 
matrix inversion. 

The inner product (/,^) is a scalar with the following properties: 

{f.g) = {gj) (3.2) 
{af + ^g,h) = a{f,h)+P{g,h) (3.3) 

here, a and /3 are scalar constants and complex conjugation is denoted by *. The 
adjoint operator L''' is defined as follows: 

{Lf,g) = {f,L^g) (3.5) 

If the operator and its adjoint are the same, then the operator is called "self-adjoint". 
Since the adjoint operator depends on the definition of the inner product and there 
is no unique way of defining the inner product; the inner product is mostly arranged 
such that the operator becomes self adjoint. A general integral expression for the inner 
product can be written as 

{f:g)= tw{x)f{x)g{x)dx (3.6) 

J a 

with the "weighting function" w(x). 

Equations I3.2H3.6I are the fundamental objects required to formulate the method of 
moments. Moving to the method of moments, let the function / of the Equation 13.11 
be expanded on the basis functions /i,/2,. . .n of the space defined by L with the 
expansions coefficients tti, tti, . . . , a„: 



f = Y.^nfn (3.7) 



n 
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Inserting this expression into 13.1 [ one gets: 



Y,CCnLfn=g 

n 

Finally, applying inner products with the weighting functions wi,W2 
both sides of 13.81 one ends up with m equations in the following form: 

n 

which can be re-written in the matrix form as: 

\l'mn\ \^n\ \Sni\ 



, . . . , Yvm 



(3.8) 



Wm to the 



(3.9) 



(3.10) 



with 



(wi,L/i) (wi,L/2) 

(w2,^/l) {w2:Lf2) 



\an\ 



[gn 



052 



(W2, 



(3.11) 



(3.12) 



(3.13) 



By the help of the matrix equation 13.101 if the expansion coefficients, a„, are known, 
it is possible to calculate the source, g, which generates that field; and also if the 
1-matrix given with 13.261 is non-singular, hence its inverse exists, one can calculate the 
expansion coefficients corresponding to a source via the inverse equation: 



[CCn] = [Ln] ^ [g, 



(3.14) 



So far, nothing has been stated about the boundaries and the weighting function of the 
inner product and the basis functions for the expansion of f and in fact, proper choice 
of these according to the nature of the problem of interest is essential in the application 
of this method. 

One approximation that is going to be used for the Poisson equation is the so-called 
point matching technique. Within the point matching approximations, one chooses the 



21 



weighting functions such that l3.1l is satisfied on a discrete set of points of interest. This 
corresponds to choosing the weighting functions simply as Dirac delta functions: 

Wm=5{x-Xm) (3.15) 

Another approximation method that is also going to be used is the choice of 
subsectional bases. In this technique, each basis function /„ is non-zero over a certain 
^th j-ggion and zero, outside that region; thus, each expansion coefficient a„ is effective 
over a certain region. In general, subsectional bases can be expressed as follows: 

r ( \ _ j f{x) if X is in the n"^ region 
/«W — I Q if X is outside the n^'' region 

These concepts of subsectional bases and point matching are going to be used in the 
solution of the Poisson equation which is explained in detail in the next subsection. 

3.1.2 Solution of the Poisson equation using moment methods 

Electrical potential in two dimensions is determined with the Poisson equation which, 
in differential form, can be written as, 

V2y = --5(z) (3.17) 
e 

here, is the Laplace operator, a is the surface charge density at z = and e is 
the electric permittivity of the medium. Dirac delta function of z trivally transforms 
the volume charge distribution to the surface charge distribution on z = 0. Assuming 
a conducting plate lying on the x-y plane of the coordinate system, the well known 
integral solution for the electrical potential on x-y plane is: 

V{x,y) = [[ y=I^^l===dx'dy' (3.18) 

J J 47te^/{x-x')^ + {y-yy 

In the operator formalism, this equation can be written as, 

Lo{x\y')=V{x,y) (3.19) 

where, 

dx'dy' 



, (3.20) 

4;re ^{x-xlY + iy-y'Y 

Equation 13.191 is the equation that is going to be written in the discrete form via the 
method of moments. A rectangular surface with a width of and length of ly, divided 
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Figure 3.1: Ix by ly rectangular surface with Hx x Uy subsections. 

in Hx subsections in x and Uy subsections in y is depicted in Figure 3.1. As it is shown 
in the Figure, each subsection is labeled with the numbers from 1 to x starting 
from the cell on the bottom left, to the cell on the top right. 

A simple approximation is to expand the surface charge distribution over the constant 
functions which are only nonzero in one cell, namely: 

f ~ (^nfn 

n 

1 if(x' ,y')isinthen^^cell 



fn{x ,y ) 



(3.21) 
(3.22) 



if (x' , y') is outside the ri^ cell 
This, physically, corresponds to assuming that the surface charge density over the area 

of the fi^ cell is constant and equal to a„. The second approximation to be done is the 

point matching approximation which is mentioned in the previous section. Denoting 

the center coordinates of each cell with (xm, jm), if the weighting function is chosen as 

Wm = 5 - Xm)b {y - y,n) (3.23) 

one can write the approximation of Equation 13.191 in the form of 13.101 the matrix 
elements Imn and the elements of the g vector can be written as follows (note that the 
function g of [XT] corresponds to j) of l3.19l) : 

Imn = {Wm.Lfn) = J J dxdyW„rLfn 

= JJ dxdyd{x-Xm)d{y-ym) 
dx'dy' 



inASn 47t£^/{x-x')^ + {y-y')^ 
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Figure 3.2: Illustration of the Imn and rifh cells and the Rmn vector. 



In 



dx'dy' 



nASn 47t£^/{x^ -X'Y + {ym -y'Y 

gm = {wm,g) = JJ dxdyd{x-Xm)d{y-yni)V{x,y) 

gm ~ ^ {p^m 5 ym ) 



(3.24) 



(3.25) 



Expression "in A5„" under the integral sign in 13.241 indicates that the integral should 
be carried out through the area of the n^^'^ cell. This is because the basis functions /„ 
13.221 are zero outside the n'^'^ cell. 



Physical interpretation of the above equations is intuitive: Imn given with the 13.241 is 
the influence of the unit surface charge distribution of the n'^ cell, on the electrical 
potential at the center of the m!^'^ cell. 

One last, harmless, approximation can be done in the calculation of the integral 
for different cells is to assume that the distance between the cells (square root in the 
denominator of l3.24l) is constant and equal to the distance between center of the cells. 
Denoting the area of the n'^ cell with A5„, matrix element can be written as, 

A5„ 



where the distance between m''^ and n^'^ cells, i?,,,,, (see Figure 3.2) is: 



th 



AneR„ 



(3.26) 



m ynj 



(3.27) 



Diagonal elements of the 1-matrix /„„ given with 13.281 this physically corresponds to 
the influence of the unit charge distribution through the ri^ surface on the electrical 
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potential at its center. Denoting the unit lengths of the cells in x and y direction with 
(^x = h/ and Uy = ly/uy, integral is computed by the help of Wolfram Mathematica 
Online Integrator ll63l as: 

1 



In 



In 



Ux In 



rax/2 i-ay/2 

dx I dy 




(3.28) 



Since a„ and have physical meanings in the approximate solution of the Poisson 
equation explained above, the approximate equation can be rewritten with the new 
symbols according to their physical meanings as follows: 

[l,nn][On] = [V,n] (3.29) 

Since the physical meaning of a„ is the value of the constant charge distribution on the 
Hth cell it is denoted by a„ instead. Similarly, since the physical meaning of the gm is 
the magnitude of the electric potential at the center of the mth cell, it is denoted by Vm- 



Combining 13.261 and I3.28[ elements of the 1-matrix can be written as: 



1 



^Jaj^af.-ay 



+ a,ln(^^±^ 



(3.30) 



m = n 



It should be noted that the elements of the 1-matrix only depends on the geometrical 
parameters, hence, once the elements of the 1-matrix are calculated, computing the 
electrical potential corresponding to a charge distribution is only multiplications and 
additions. This is desired, as through the transient simulation, this calculation has to 
be done many times. 

It is also possible to calculate the charge distribution corresponding to a certain 
potential by using the inverse equation of l3.29t 

[On] = [l,nnr'[Vm] (3.31) 

For a constant potential, charge distribution and, hence, the total charge can also be 
calculated using 13.311 One can find the capacitance of the single plate by calculating 
the total charge corresponding to the unit potential. 
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3.2 Chani 

The developed simulation tool for the charge transport and discharge on a rectangular 
surface is named "Chani". This is due to both words "charge" and "Chani" sharing the 
first two letters "ch" and referral to the name of a character from the science fiction 
novel Dune li64ll . 

Chani is developed to be run on the ROOT Framework f65] with which running C++ 
code as a macro (without compiling) is possible via CINT, C++ interpreter If66l . Chani 
also takes the advantage of the ROOT vector, matrix and histogram classes and the 
libraries for visualization. 

The main functions of the Chani are declared in the header file, Chani.h. These main 
functions and their jobs are explained as follows: 

parameters(lx, nx, ly, ny, Nm, sigmas): Defines the surface with dimensions Ix and ly 
with surface resistivity sigmas. Tells the simulator to divide the surface in x and y to 
nx and ny respectively. States that the 1 matrix is going to be calculated in Nm x Nm 
parts. 

computelmn(): Calculates the 1-matrix. 

addConnector(x, y, R): Adds a ground connection with connector resistance R to the 
nearest cell to (x, y). 

addCharge(x, y, q, tStep): Adds q amount of charge to the nearest cell to (x, y) at the 
time instance tStep. 

transient l(tStepi, nofSteps, Deltat): Performs the transient calculation starting from 
the time step "tStepi" for "nofSteps" steps of "Deltat" seconds. Electrical potential 
and charge distribution objects are written to and read from hard drive at each step. 
Suitable for large-scale calculations, where, memory is not enough to hold all of the 
objects. 

transient2(tStepi, nofSteps, Deltat): Performs the transient calculation starting from 
the time step "tStepi" for "nofSteps" steps of "Deltat" seconds. Electrical potential 
and charge distribution objects are held in the memory. Suitable for small-scale 
calculations. 

capacitancei): Calculates and outputs the capacitance of the surface. 
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Figure 3.3: Resistances and the currents between neighboring cells. 

getV( i): Gets the i-th element of electrical potential vector in function transient 1(). 
getQtotal(tStep): Gets the total amount of charge at the time step tStep. 
getQxy(x, y, tStep): Gets the amount of charge on the nearest cell to (x, y) at the time 
instance tStep. 

Header file, main functions and main classes are given in Appendix A. 

A rectangular surface is handled by the simulator as explained in the previous section 
(see Figure 3.1). Through the transient simulation electric potential of each cell and 
the currents between neighboring cells are calculated. Surface resistivity (Ps), assumed 
to be constant over the surface, is an input of the simulator. For a cell with the label /, 
equivalent resistance between the i'^ cell and the cell on its right (cell number: /+ 1) 
is: 

Rx = {ps X ax) /ay (3.32) 
Similarly, the resistance between f'^ cell and the upper cell (cell number: % + /) is: 

Ry = {ps X ay)/ax (3.33) 

Hence the currents between these cells are simply calculated from Ohm's law as: 

Ix, = {Vi-Vi+i)/Rx (3.34) 

Iy, = {Vi-Vn,+i)/Ry (3.35) 

Here, V/'s are electrical potential. These currents and resistances are illustrated in 
Figure 3.3. 

These currents are assumed to be constant during one time step. Hence the amount 
of the migrating charge (A^) from on cell to the other cell is simply calculated by 
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multiplying the current with the length of the time step (At): 

Aq = IxAt (3.36) 

Finally, resistances at the ground connection points and the cells that the ground 
connection is made, are: 

Rg = Rc+ [ps X (a^(y)/2)]/a^(^) (3.37) 

If the connector is in x (y) direction. In the equation above, Rc is the "connector 
resistance" which is the resistance between the power supply and the connection point, 
and it is summed with the resistance coming from the half-length of the connection 
cell. Discharge current and the amount of charge removed is calculated similarly. 

In the next section, the charge transport simulation on a 2 x 2cm^ surface carried out 
by Chani is explained. 



2 cm 



J 

•- ^ 

2 cm 

Figure 3.4: 2 x 2cm^ surface with total number of 225 cells. 
3.3 An Example 

Illustration of the surface on which the charge transport is simulated using Chani is in 
Figure 3.4. It is a 2 x 2cm^ square surface with the constant resistivity of lO^Q/D. 
In the simulator, it is divided to 15 in both directions hence the total number of cells 
for the simulation is 225. The surface is grounded from the left end of it without any 
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connector resistance and initially 10 elementary charges arrive on the center and the 
neighboring eight cells of the surface. 

Transient simulation with Chani is performed for 2000 time steps with one step being 
10^*^^ long, thus physically for 200ns. Surface definition and simulation files are 
given in Appendix B & B.l respectively. 

Change in the total amount of charge with time is plotted in Figure 3.5 where 
exponential-like discharge is clearly seen. It should be noted that no capacitance 
parameter is input to the simulation but the effect of self-capacitance of the plate 
intrinsically calculated and thus RC-like discharge behavior is obtained. 



Total charge versus time 

XI 0-'^ 




Time (s) 

Figure 3.5: Total charge on the surface versus time. 

Figure 3.6 shows the surface charge density at four time instances. It can be seen 
that the charge, initially at the center is spread very quickly and at t = 20 ns, charge 
distribution becomes almost uniform except the edges and the ground connection point 
as expected. 

3.4 Self- Consistency Tests 

Two self-consistency tests are performed in order to proove the convergence of the 
calculations and to estimate the order of errors. In this section, these tests are explained 
and example calculation results on the surface studied in the previous section is 
presented. 
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Surface charge density at t = ns Surface charge density at t = 1 ns 




Surface charge density at t = 1 ns Surface charge density at t = 20 ns 




Figure 3.6: Color plot of the surface charge density at four different times. 
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3.4.1 Capacitance versus number of cells 

One way to check the accuracy of the calculations is to check the capacitance of the 
surface with the increasing number of cells. One should expect to find capacitance 
value converging to a certain real value as the number of divisions increases. Although 
it is not directly related with the transient calculations, convergence of capacitance is a 
reference point to check whether the simulation is meaningful. 

As an example, capacitance of the 2 x Icm^ surface is iteratively calculated while 
increasing the number of divisions. After each iteration relative error on capacitance 
is calculated as follows: 

Error = (3.38) 

Cprevious 

Figure [377] shows the variation of capacitance against number of divisions used in the 
calculation and the variation of error defined with Equation l3.38l is shown in Figure lXSl 
Iterative calculation of capacitance is performed until the relative error goes below the 
tolerance 0.001. Source file can be found in Appendix B.2. The tolerance is satisfied 
with 17 divisions in each direction hence 289 total number of cells. It is clearly seen 
from the figures 3.7 and 3.8 that the capacitance approaches to a certain value with 
decreasing amount of error as the number of cells increases. 

Capacitance vs Number of divisions 
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Figure 3.7: Capacitance versus number of cells. 
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Relative Error vs Number of divisions 
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Figure 3.8: Relative error versus number of cells. 



Another "sanity check" may be to check the charge distribution corresponding to 
constant potential. capacitance() function of Chani also saves the charge distribution 
corresponding to the unit potential. Surface charge density distribution of the 2 x 2cm^ 
corresponding to unit consant potential is shown in Figure [X9l It is clearly seen that the 
more charge stands on the comers than the central points and this is indeed expected 
since the charge on a conductor accumulates on the sharp edges. 



Surface charge distribution corresponding to unit potential 




Vol -0.008 -0.006 -0.004 -0.002 
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Figure 3.9: Surface charge density distribution corresponding to the unit constant 
potential. 



3.4.2 Discharge time versus number of cells 

As a second self-consistency check, the time needed for 90% of the inital charge to 
be removed is calculated iteratively with the increasing number of subdivisions. The 
number of cells is increased until the relative error between the previous discharge 
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time is less than 0.001. Denoting the discharge time with Relative error is defined 
as follows: 

Error ^d, previous iteration I (3 39) 

^d, previous iteration 

2 X 2cm^ surface with lO^fl/D surface resistivity is first divided by 7 in each direction 
resulting total 149 cells and initially 10^ elementary charge is added in the center 
similar to the simulation described previously (see Figure 3.2). The surface is grounded 
from the left end again in the same way as shown in Figure 3.2. Transient simulation 
is performed until the 90% of the initial charge is removed and then the number of 
divisions in each direction is increased by 2 and all other steps are repeated. Related 
code is in Appendix B.3. This process is iterated until the relative error (as defined 
with Equation 3.39) between discharge times from two adjacent iterations becomes 
less than 0.001. Why the number of divisions are increased by 2 in every iteration is 
that if they were increased by 1 , number of divisions were going to be even in every 
2 iterations and that is not desired because for even number of divisions there is no 
"central cell", thus, it is not possible to put the initial charge at the same coordinates 
for even and odd cases. 

Discharge time vs Number of divisions 
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Figure 3.10: Discharge time versus number of cells. 

The 0.001 tolerance on the relative error is satisfied for 21 divisions in each direction 
(total 441 cells). Variation of discharge time and the relative error with the increasing 
number of cells is plotted on Figure [37TOl and Figure [37TT] respectively. It can be clearly 
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seen that discharge time asymptotically approaches to a value with decreasing amount 
of relative error as the number of divisions increases. 



Relative Error vs Number of divisions 
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Figure 3.11: Relative error in discharge time versus number of cells. 



3.5 An Analytically Solvable Case 

Relaxation of a periodic charge distribution on an infinite conducting plane is described 
by a simple analytical solution. In the next subsection, this analytical solution is 
explained; and in the following subsection, results of several simulations performed 
by Chani, where a conducting square plate had initially a periodic charge distribution, 
are presented. 

3.5.1 Infinite conducting plane with periodic charge distribution 

The periodic surface charge distribution on an infinite plane lying on z = is described 
with the following equation: 

o = Oocos(kx)exp{-t/t) (3.40) 

Here, periodicity is determined by the cosine function and the exponential is the 
relaxation term. Boundary condition on z = from Gauss' Law reads: 

E+-E- = ^ (3.41) 

Co 

dz dz Co 
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In these equations, lower index z denotes the z-component of the electrical field and 
upper indices + and — denotes "just above" and "just below" the z = surface 
respectively. Electrical potential satisfying Equation 13.421 for the charge distribution 
given with 13.401 is: 

V = -^cos(fc;c)exp(-A:|z|)exp(-?/T) (3.43) 

Inserting this to l3.42l one can see that the boundary condition at z = is indeed satisfied: 

k^^^ cos{kx) exp{—t / t) — (— fc)^^^cos(fcjc)exp(— ?/t) = ■^cos(fcx)exp(— ?/t) 

a 

Continuity equation in two dimensions can be written as, 

d(5 

— + V-K = (3.44) 

ot 

Here, K denotes the surface current density and V contains derivatives in x and y. 
Ohm's law for surface currents is: 

K = — E (3.45) 

Ps 

Here is the surface resistivity. Writing Ohm's law in terms of electrical potential 
and substituting it into Equation 13.441 one has, 

(9(7 1 9 

VV = (3.46) 

dt p. 

Finally, inserting the surface charge distribution and electrical potential at z = given 
with the Equations 13.401 and I3.43l respectively into Equation 13.461 one has, 

— -aocos(fcjc)exp(— ?/t) +^^r^ — cos(fcx)exp(— ?/t) = 



which yields. 



(3.47) 

k 



This time constant T characterizes the relaxation time of the periodic charge 
distribution. In the next subsection, x is extracted from the fits to the simulation results 
from Chani and compared with the values given with 13.47] 
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3.5.2 Simulations of the charge relaxation 

In the simulations of this section (see Appendix C for the source files), a square plane 
with no ground connection, initially has a periodic charge distribution given with 13.401 
Transient calculation is performed and at each subcell, variation of the amount of 
charge with time is obtained. Decay of the charge in the middle is fitted with an 
exponential and the corresponding time constant is extracted. As an example, decay of 
the charge in the middle of the surface from the simulation of 2 x 2 cm^ is plotted in 
Figure [3TT2I Color plots of the surface charge distribution at four time instances from 
the same simulation is shown in Figure [BTSI 

Charge at (0,0) versus time 




Time (ns) 



Figure 3.12: Decay of the total charge in the center of a square plate which initially 
had a periodic charge distribution. Red curve is the exponential fit. 

Time constants corresponding to the relaxation of periodically distributed charge 
are obtained from the simulations performed by Chani for squares with different 
sidelengths. In these simulations, the wavenumber corresponding to the charge 
distribution given with 13.401 is fixed to lOcm^^ and the surface resistivity is lO^tl/D. 
The time step is 0.01 ns. Results are presented in Table [BTI where / is the sidelength 
of the square plate, n is the number of divisions in both sides (total number of subcells 
is nx n) and Tchani is the time constant that is obtained from the exponential fit to the 
results of the simulation (as shown in Figure [3TT21) . As it can be seen from Figure l3T3l 
periodic distribution of the charge on the surface is lost with increasing time which is 
not the case when the plane is infinite. For this reason, exponential is fit to only the 
first 2 ns of the simulation assuming that during this time interval the behavior should 
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be similar to the infinite plane solution. For k = lOcm ^ and = \O^Q./D, the time 
constant for the infinite plane is calculated from|3]47]is Tinfinite plane = 

1.7708/15. The 

%Error is defined as follows: 

„ „ I ^infinite plane " - T'Chani , ,^ .^^ 

%Error = x 100 (3.48) 

^infinite plane 

As it can be seen from the Table 13.11 error between the time constants from the 
theory and Chani is less than 3.5% for 5 different lengths. One might expect to see a 
monotonic decrease of the error as the size increases assuming that the behavior would 
become much similar to that of the infinite plane. However, this is not the case with 
these results indicating that there must be some other effect related with the simulation 
conditions. 

Table 3.1: Relaxation time constants for different sidelengths. 



l(cm) 


n 


TChani (ns) 


% Error 


1 


25 


1.7363 


2.0 


2 


49 


1.8146 


2.5 


3 


75 


1.8320 


3.5 


4 


99 


1.8300 


3.3 


5 


125 


1.8224 


2.9 



The same simulation procedure is followed, this time, keeping the sidelength constant 
and equal to 5 cm while varying the number of divisions (n). Results obtained from 
the fit and errors are listed in Table 13.21 The wave number, k, is again fixed to 1000 
and the time step is 0.01 ns. 

Table 3.2: Relaxation time constants of the square plate with 5 cm sidelength with 
different number of divisions. 



n 


'TChani (ns) 


%Error 


9 


14.7735 


734 


13 


4.6277 


161 


17 


3.4186 


93 


25 


2.5793 


46 


33 


2.2641 


28 


41 


2.1092 


19 


49 


2.0205 


14 


57 


1.9643 


11 


65 


1.9262 


9 
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Surface charge density at t = ns 



Surface charge density at t = 1 ns 
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Surface charge density at t = 7 ns 
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Figure 3.13: Color plot of the surface charge density at four different times. 
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Results in Table 13.21 shows that the error with the theory decreases with increasing 
number of divisions as expected. In the first two rows of the Table 13.21 the error is 
above 100%. In fact, these two calculations do not make sense at all because they do 
not satisfy the "Nyquist rate" Il67l . Initial charge of the surface was periodic with the 
wave number k = lOOOm^ which corresponds to the spatial frequency: 

f = k/2n=l59.l5m-^ (3.49) 

In essence, Nyquist - Shannon sampling theorem tells that a signal with a bandwith 
B, should be sampled with a sampling frequency greater than 2B in order to be 
able to reconstruct the original signal without aliasing. For the simulations of this 
section, charge distribution function [3.401 and the subcells of Chani can be thought as 
the original signal and the sampling points respectively. Sampling frequency can be 
expressed as: 

f^ = {l/n)-^ =n/l (3.50) 
Where 1 is sidelength and n is number of divisions. Nyquist criteria reads, 

/, > 2/ (3.51) 
n>2x 159.15m"^x 0.05m =15.9 (3.52) 

This simple calculation tells that if the number of divisions is chosen less than 16, 
initial charge distribution is not properly understood by the simulator hence result it 
gives is unreasonable. 

In the Table [X2l 93% error is given for n = 11 which is just above the Nyquist limit 
indicating that a higher sampling frequency is required. Error with the theory falls 
below 10% when n = 65 for which the sampling frequency is more than 4 times the 
Nyquist frequency. 

3.6 Charge Transport on a Resistive Strip 

Up to this point, presented calculations were on a square shaped surface. However, 
main motivation of this project was to develop a tool for the optimization of the 
resistive structures in the resistive-anode micromegas detectors being developed by 
MAMMA group and these resistive strips are rectangular. For example, prototypes 
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presented in Q have 10cm by 150/im resistive strips with surface resistivities 
0.030,0.075, and0.0075fl/n. As an example, a transient simulation of charge 
transport and discharge is performed with the surface parameters from Q. 

Proper choice of cell dimensions becomes essentially for calculations on rectangular 
surfaces with big length-to-width ratios. Several trials showed that divisions must be 
made such that the unit cells of the simulations have a square-like shape. Unreasonable 
results for charge distribution is obtained when rectangular unit cells with high 
length-to- width ratios are chosen. 



r 



ly = 1 50 ptm 



Hy = 1 



Ix = 10 cm, Hx = 667 
Figure 3.14: Illustration of resistive strips with divisions. 



10cm by 150/im rectangular surface, illustrated on Figure [3T4l is divided by 667 in 
X and not divided in y. Initially 10"* elementary charge is put in the middle and the 
surface is grounded from the left end as it is shown in Figure [3T4l Surface resistivity 
is set to 0.03 n/D and transient calculation is performed for 2,000,000 steps with 0.1 
ns time steps. See Appendix D for source files. 
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Figure 3.15: Discharge plot for the rectangular surface. 
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Change in the amount of total charge with time is shown on Figure [3TT51 As it is shown 
in the figure, the time needed for the removal of most of the initial charge is 30 {is. The 
frequency corresponding to this period is l/SOjAs ^ 33 kHz. Thus one should expect 
a charge up on the resistive strips when the detector is put under an incoming particle 
rate greater than this frequency. 

Color plots of the surface charge density at four simulation instances are given in 
Figure l3TT6l Spread of the initial charge can be seen on the plot on top right. 

3.7 A Large Scale Calculation 

Last simulation that is going to be presented in this section is again a resistive structure 
from one of the prototypes (R14) build by MAMMA group. Similar to the previous 
simulations, spread of the initial charge in the center is studied. This surface is a 
100 X 18 mm^ rectangular surface with 0.24 MQ./n surface resistivity. The different 
thing about this simulation is it has a large number of cells. It is divided by 51 in one 
direction and 279 in other resulting 14229 cells. 

Calculation with 14229 cell corresponds to a 14229 by 14229 1-matrix and this 
is a huge matrix dimension to hold in the memory. Chani does the calculations 
with big 1-matrices by dividing the matrix to submatrices. After that, calculation 
of the potentials are done via sub-vectors, written in the harddrive and read for the 
current calculations, this is done by the transient 1() function which works slower than 
transient2( ) but has higher memory capabilities. Number of submatrices is a user input 
parameter and should be carefully choosen. 

Transient simulation on the 100 x 18 mm^ surface is performed for 600 steps with 
time step 0.1 ns. Simulation files are given in Appendix E. Color plot of the surface 
charge density is given on Figure [JTTl where charge spread can be seen. 
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Figure 3.17: Color plot of the surface charge density at four different times. 

m 
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4. CONCLUSIONS 



Main goal of the study that is presented in this thesis was to understand the charge 
spread and discharge dynamics of the electrical charge on a rectangular surface of high 
resistivity in order to help the design of the micromegas detectors with resistive anodes. 
For this purpose, a simulator is developed and applied on several cases with different 
dimensions and resistivities. 

It is shown via self -consistency tests in sections 3.4.1 and 3.4.2 that the simulations are 
convergent if the subdivision parameters are properly chosen. 

Relaxation of periodic charge on square shaped surfaces is studied and the simulation 
results are compared with the theoretical expectations in section 3.5. With properly 
chosen parameters, error between values of the relaxation time constants from theory 
and Chani were below 3.5%. 

The calculation on a resistive strip which is presented in the Section 3.6 demonstrates 
that the starting motivation, calculating the charge transport on a resistive strip in 
micromegas chambers being developed by MAMMA group is indeed achieved. This 
calculation also provides the fact that the simulator works for very long time steps 
(2,000,000) without any memory overflow or instability. 

In Section 3.7, a simulation with 14229 cells is presented, proving that the simulator 
can handle very large matrix operations. 

In conclusion, a working simulator for the charge transport on rectangular surfaces is 
developed. It is possible to calculate the time needed for the total discharge and the 
spread of the charge over to surface using Chani. Although it is developed to help 
the design of the gaseous ionization detectors, Chani's features can find applications 
in many fields where the understanding of charge transport properties is important. 
Specifically for particle detectors with resistive anodes, Chani can be used to guess 
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the rate capabilities of the detectors from the discharge times and possible cross-talks 
between the readout strips due to the charge spread. 
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APPENDIX A Chani.h header file, main function and classes 

#ifndef _CHANIH__ 
# define _CHANIH__ 

/^Header file of Chani (Charge transport simulator for rectangular s 
urfaces ) . 

* su rfa ces 

* Created: Apr 25, 2012 by Burak Budanur ( burakbudanur® gmail . c 
om) 

* Last edited: May 21, 2012 by Burak Budanur ( burakbudanur© gma 
i I . com ) 

*/ 

#include <iostream> 

//Root libs: 

#include "TMath.h" 

#include "TMatrixD.h" 

#include "TVectorD.h" 

#include <cmath> 

#include "TFile.h" 

#include "THl.h" 

#include "TH2.h" 

//Chani functions : 

#include " functions / parameters .C" 

#include " functions / computelmn .C" 

#include "functions / addConnector .C" 

#include "functions /addCharge .C" 

#include " functions /transientl .C" 

#include " functions / transient2 .C" 

#include " functions /getV .C" 

#include "functions / getQtotal .C" 

#include " functions / capacitance .C" 

#include " functions /getQxy .C" 

//Chani classes : 

#iiiclude " classes /Readp .C" 

#include " classes / Findij .C" 

//Global constants 

#include " globalConsts . h" 

using std : : cout ; 
using std : : endl ; 

void parameters ( double Ix , int nx , double ly , int ny , int Nm, double 

sigmas ) ; 

void computelmnO ; 

void addConnector ( double x, double y, double R) ; 

void addCharge ( double x, double y, double q, int tStep); 
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void transientl ( int tStepi , int nofSteps , double Deltat); 

void transient2 ( int tStepi , int nofSteps, double Deltat); 

double capacitance () ; 

double getV ( int i ) ; 

double getQtotal ( int tStep); 

double getQxy ( double x, double y, int tStep); 

#endif // __CHANIH__ 

#ifndef _GLOBALCONSTSH_ 
# define _GLOBALCONSTSH_ 

#include "TMath.h" 

const double pi = TMath : : Pi ( ) ; 

const double epsilonO = 8.8541878176e-12; // F / m 
const char * dir = " /$ROOTSYS/ burak / ChaniOv3sine / output /" ; / 
/ DIRECTORY TO HOLD OUTPUT FILES . EDIT THIS LINE ACCORDING T 
O YOUR WORKING DIRECTORY . 

#endif 

#include " . ./ Chani .h" 
#include " . . / globalConsts . h" 

void parameters ( double Ix , int nx , double ly , int ny , int Nm, double 

sigmas) { 



TString fnameg = TString :: Format ( "%sparameters . root " , dir); 

TString fnamev = TString :: Format ( "%sxy . root " , dir); 

TString fnameh = TString :: Format ( "%sh . root " , dir); 

double ax = Ix/nx; // unit length in x 
double ay = ly/ny; // unit length in y 

int N = nx*ny ; 

std::cout « "The total number of subsections is N = " « N « std : 
: endl ; 

std::cout « "The 1— matrix is going to be calculated in "; 
std::cout « Nm*Nm « " parts." « std :: endl ; 



// Assigning the x & y coordinates to the subsections 

TVectorD x(l ,N) ; 
TVectorD y(l ,N) ; 

if (nx%2 == && ny%2 == 0) { 

for (int j = 1; j < ny + 1; j++){ 
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for (int i=l; i < nx + 1; i++){ 

x( i+nx*(j — 1)) = — ax*(nx/2) + ax/2 + (i— l)*ax; 
y( i+nx*(j — 1)) = — ay*(ny/2) + ay/2 + (j— l)*ay; 

} 

} 

} 

else if (nx%2 == 1 && ny%2 == 1){ 
for (int j = 1; j < ny + 1; j ++) { 
for (int i = l; i < nx + 1; i++){ 

X ( i +nx * ( j — 1) ) = — ax * (( nx — 1) /2) + (i— l)*ax; 
y( i+nx*(j -1)) = -ay*((ny-l)/2) + (j-l)*ay; 

} 

} 

} 

else { 

std : : cout « "number of subsections for x and y should be both ev 
en or both odd." « std : : endl ; 
return ; 

} 

TH2 * sigmahO ; 

sigmahO = new TH2D( " sigmaO (x ,y ) " , "Surface charge density at t = 0" 
, nx, -lx/2, lx/2, ny, -ly/2, ly/2) ; 

TH2 *connh; 

connh = new TH2D( " conn (x , y ) " , "Connectors histogram", nx , —lx/2, Ix 
12, ny, -ly/2, ly/2) ; 

for (int j = 1; j < ny+1; j ++) { 
for (int i = l; i < nx+1; i++){ 

connh— >S etCellContent ( i , j, le20); // Assign a nonzero dummy valu 
e to the histogram contents 

} 

} 



TFile f (fnameg ," recreate ") ; // open the output file. 

THl *params; // The histogram to hold Ix , nx , ax, ly , ny , ay, Nm va 
lues in 

params = new TH1D( " params " , " Ix , nx , ax, ly , ny , ay, Nm",8,l,8); 



params — >SetBinContent ( 1 , Ix ) 

params —>SetBinContent (2 , nx) 

params —>SetBinContent (3 , ax) 

params —>SetBinContent (4 , ly) 

params —>SetBinContent (5 , ny) 

params —>SetBinContent (6 , ay) 

params —>SetBinContent (7 , Nm) 
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params— >SetBinContent (8 , sigmas) ; 

params— >Write ; // Write params histogram in file 
delete params ; 

f. Close 0; // Close the output file 

TFile f2 ( fnamev ," recreate ") ; // open the output file for x & y. 

X . Write ( "x— vector ") ; // Write x coordinates in the file 
y . Write ( "y— vector ") ; // Write y coordinates in the file 

f2 . Close ; 

TFile f 3 ( fnameh ," recreate ") ; // open the output file for histograms 

sigmahO — >Write ( " sigmahO " ) ; // Blank histogram of the surface 
connh— >Write ( " connh" ) ; // Connector histogram of the surface 

f3 . Close ; 

delete connh; 
delete sigmahO; 

} 

#include " . . / Chani .h" 
#include " . . / globalConsts . h" 
#include "../ classes /Readp . h" 

void computelmn { 

Readp rp ; 

double Ix = rp.GetlxO; 
double ly = rp.GetlyO; 
double ax = rp . Getax () ; 
double ay = rp.GetayO; 
int nx = rp.GetnxO; 
int ny = rp . Getny () ; 
int Nm = rp .GetNm() ; 

int N = nx*ny ; 

TString fnamev = TString :: Format ( "%sxy . root " , dir); 

// Reading of x & y vectors: 
TFile f ( fnamev ) ; 

TVectorD *xp; // Pointer for x— vector 
f . GetObjectC'x-vector" ,xp) ; 
TVectorD x(l ,N) ; 
X = *xp; 
delete xp; 

TVectorD *yp; // Pointer for y— vector 
f . GetObjectC'y-vector" ,yp) ; 
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TVectorD y(l ,N) ; 

y = *yp; 
delete yp ; 

f. Close ; 

double DeltaSn = ax*ay; //Area of a subsection 

std::cout « "Area of a subsection DeltaSn = " « DeltaSn « 

std : : endl ; 

std : : cout « "Calculation of the l(m,n) matrix ..." « std :: endl ; 

int nsub = N / Nm; 

TMatrixD Isub ( 1 , nsub , 1 , nsub ) ; 

for (int i = 1 ; i <= Nm; i ++) { 
for (int j = 1; j <= Nm; j++){ 

TString fname = TString :: Format ( "%slmatrix%i . root ", dir , i+(j-l)*N 

m); 

TFile f (fname ," recreate ") ; // open the output file. 

for (int m = (j — l)*nsub +1 ; m <= j*nsub; m++){ 
for (int n = (i — l)*nsub + l; n<= i*nsub; n++){ 

int mm = m — (j — l)*nsub ; 
int nm = n — (i— l)*nsub; 

if (n == m){ 

Isub (mm, nm) = (l/(4*pi*epsilon0)) 

*(ax*log (( sqrt (pow(ay ,2) + pow(ax,2)) + ay ) /( sqrt (pow(ay ,2) 
+ pow (ax, 2)) — ay)) 

+ ay*log (( sqrt (pow(ax ,2) + pow(ay,2)) + ax) /( sqrt (pow(ax ,2) 
+ pow(ay ,2) ) - ax) ) ) ; 
} 

else { 

lsub(mm,nm) = DeltaSn /(4* pi* epsilonO* sqrt (pow((x(m)—x(n) ) ,2)+p 
ow((y(m)-y(n)) ,2))) ; 
} 

// std :: cout « « m « ", " « n « ") = " « Isub (nm,nm) 

« " " ; 
} 

} 

char label [12]; 

sprintf (label , " l-matrix%i " , i +(j -l)*^&n) ; 
Isub . Write (label ) ; 

f. Close 0; // Close the output file 

} 

} 
} 

#include " . ./ Chani .h" 
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#include " . . / globalConsts . h" 
#include "../ classes / Findij . h" 

void addConnector ( double x, double y, double R) { 

//Names of the files to be read 

TString fnameh = TString :: Format ( "%sh . root " , dir); 

Findij f ij (x,y) ; 
int i = f ij . Geti () ; 
int j = fij .Getj ; 



TFile f (fnameh, "UPDATE"); // open the output file for histograms 

TH2 *connh; // The histogram to hold Ix , nx , ax, ly , ny , ay values 
in 

f . GetObject("connh" ,connh) ; 

connh— >SetCellContent(i , j , R) ; // Write the resistance value in th 
e connector cell 
connh— >Write (" connh" ) ; // Connector histogram of the surface 

delete connh; 

f .CloseO ; 

} 

#include " . ./ Chani .h" 
#include "../ globalConsts . h" 
#include "../ clas ses / Findij . h" 

void addCharge ( double x, double y, double q, int tStep){ 

// File name : 

TString fnameh = TString :: Format ( "%sh . root " , dir); 

Findij fij (x ,y) ; 
int i = f ij . Geti () ; 
int j = fij .Getj () ; 
int nx = fij .Getnx() ; 
int ny = fij .Getny() ; 
double ax = fij .Getax(); 
double ay = fij.GetayO; 

double sigmaadd = q/(ax*ay); // additional surface charge density 

TFile f (fnameh, "UPDATE"); // open the output file for histograms 
// f . GetListOfKeys ( )—>P rint ( ) ; // Print to see if the objects are re 
trieved properly 

TH2 *sigmah; // The temporary surface charge density histogram 
char objName[16]; 

sprintf (objName , " sigmah%i " , tStep); 
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f . GetObject(objName , sigmah); 

sigmah— >SetCellContent( i , j , sigmah— >GetCellContent(i , j ) + sigmaadd) ; 
sigmah— >Write ( objName) ; // Write the histogram back 
f. Close ; 

} 

#include " . ./ Chani .h" 
#include "../ classes /Readp . h" 
#include "getV.C" 

void transientl ( int tStepi , int nofSteps , double Deltat){ 
Readp rp ; 

double Ix = rp.GetlxO; 

double ly = rp.GetlyO; 

double ax = rp.Getax(); 

double ay = rp . Getay () ; 

int nx = rp.GetnxO; 

int ny = rp.Getny(); 

int Nm = rp .GetNm() ; 

double sigmas = rp . Getsigmas () ; 

int N = nx*ny; // number of sub — c ells 

int nsub = N / Nm; // dimension of the sub matrices 

int Nl = Nm * Nm; // number of I— matrices 

double DeltaSn = ax*ay; // unit area 

//Create and write Vi & sigmai vectors: 

for (int i = 1; i <= Nm; i++){ 

TVectorD V(l ,nsub) ; 
TVectorD sigma ( 1 , nsub ) ; 

TString fname = TString :: Format ( "%sVsigma%i " , dir , i); 

TFile f (fname, "recreate"); // open the output file 

char labelV[5]; 
char labelsigma [ 1 0] ; 

sprintf (labelV , "V%i " , i); 
V. Write ( labelV ) ; 

sprintf ( labelsigma , " sigma%i " , i); 
sigma . Write (labelsigma) ; 

f. Close ; 

} 

//Read the initial surface charge density and connector histogram : 
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TH2 *sigmaht; // Surface charge density histogram 
sigmaht = new TH2D( " SigmahtStep (x , y ) " , "Surface charge density" , 
X, -lx/2, lx/2, ny, -ly/2, ly/2) ; 

TH2 *connh; //The connector histogram : 

connh = new TH2D( " Connectors " , "Connector histogram" , nx , —lx/2, 
x/2, ny, -ly/2, ly/2) ; 

TString fnameh = TString :: Format ( "%sh . root " , dir); 

TFile f (fnameh); // open the output file for histograms 

char objName[16]; 

sprintf (objName , " sigmah%i " , tStepi); 

TH2 * sigmahtd ; 
TH2 *connhd; 

f . GetObject (objName , sigmahtd); // Read the initial surface charge 
density histogram 
f . GetObject (" connh" , connhd) ; // Read the connector histogram 

for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 
sigmaht->SetCellContent(i , j, sigmahtd ->GetCellContent ( i , j )) ; 

} 

} 

for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 
connh— >SetCellContent ( i , j, connhd— >GetCellContent ( i , j )) ; 

} 

} 

delete sigmahtd ; 
delete connhd; 

f. Close ; 

//Write the initial surface charge density histogram to the sigma 
ectors : 

int dunomy; 

for (int isigma = 1; isigma <= Nm; isigma++){ 

TString fname = TString :: Format ( "%sVsigma%i " , dir, isigma); 

TFile f2 (fname, "UPDATE"); // open the output file 
///. GetListOfKeys ( )->P rint ( ) ; 

char labelsigma [ 1 0] ; 

sprintf ( labelsigma , " sigma%i " , isigma); 

TVectorD * sigmav ; 

f2 . GetObject( labelsigma , sigmav ) ; 
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for (int j = 1; j <= ny ; j++){ 
for (int i = 1; i <= nx; i++){ 

dummy = ((j — l)*nx+i — l)/nsub + 1; 

if (dummy == isigma){ 

// std : : cout « "cell (" « i « ", " « j « ") belongs to the 
vector "« isigma « std : : endl ; 

(* sigmav ) (( j — l)*nx+i — ( isigma — 1)* nsub ) = sigmaht— >GetCellCont 
ent(i, j); 

} 

} 

} 

sigmav —> Write (labelsigma) ; 
delete sigmav; 
f2 . Close ; 

} 

double Kt , It, Rcon ; //Surface current density, Discharge current 
TVectorD Vi(l,iisub); // i-th V vector 
TVectorD sigmavj ( 1 , nsub ) ; // j— th sigma vector 
TMatrixD Ij ( 1 , nsub , 1 , nsub ) ; // [ ( i - l)*Nnw-j ]- th 1-matrix 

// Transient loop : 

for (int it = tStepi + 1; it <= tStepi + nofSteps ; it++){ 
std :: cout « "Time step " « it « std :: endl ; 
//Calculation of the new potential: 
for (int i = 1; i <=Nm; i++){ 
Vi.ZeroO ; 

for (int j = 1; j <=Mn; j++){ 

//std:: cout « "i , j = " « i « ", " « j « std :: endl ; 

TString fname = TString :: Format ( "%sVsigma%i " , dir , j); 

TFile fl (fname); // open the output file of vectors 

TVectorD *sigmavp; // pointer for the sigma— vector 

char labelsigma [ 1 0] ; 
sprintf ( labelsigma , " sigma%i " , j); 
fl . GetObject ( labelsigma , sigmavp) ; 
sigmavj = * sigmavp; 

delete sigmavp ; 
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fl . Close ; 



/ / sigmavp — >~sigmavp ( ) ; 
// delete sigmavp ; 

TString fnamel = TString :: Format ( "%slmatrix%i . root ", dir , (i-l)*N 
mfj ) ; 

TMatrixD * Ij p ; // pointer for the (i— l)Nmfjth 1— matrix 
TFile f (fnamel); // open the output file of l(i— l)Nm+j 
//f .GetListOfKeys()->Print ; // 
char label[10]; 

sprintf ( label , " l-matrix%i " , ( i -l)*N&n+j ) ; 
f . GetObject(label , Ijp ) ; 
//*/ 
f .CloseO ; 

Ij = *ijp ; 
delete Ijp ; 

Vi = Vi + Ij * sigmavj ; 

} 

//std::cout « " Vi(38) = " « Vi(38) « std : : endl ; 

TString fname = TString :: Format ( "%sVsigma%i " , dir, i); 

TFile f (fname, "UPDATE"); // open the output file 

char labelV[5]; 

sprintf (labelV , "Woi " , i); 

Vi. Write (labelV) ; 

f .CloseO ; 



//Update the surface charge density histogram due to the currents: 
// Handling the currents in x: 

for (int j = 1; j <=ny ; j ++) { 
for (int i = 1; i<=nx — 1; i ++) { 

Kt = (( — 1) *(getV ( i+l+nx *(j — 1)) — getV ( i+nx*(j — 1)) ) * sigmas ) / ax ; 
sigmaht->SetCellContent(i , j, sigmaht->GetCellContent ( i , j ) - (Kt 
*Deltat*ay)/DeltaSn) ; 

sigmaht->SetCellContent ( i +1 , j, sigmaht->GetCellContent ( i +1 , j ) + 
(Kt*Deltat*ay)/DeltaSn) ; 

} 

} 



64 



// Handling the currents in y: 
for (int j = 1; j<=ny-l; j ++) { 
for (int i = 1; i<=nx; i++){ 

Kt = (( — 1 ) * ( getV ( i +nx* j ) — getV ( i+nx * (j — 1) ) ) * sigmas ) / ay ; 
sigmaht->SetCellContent(i , j, sigmaht->GetCellContent(i , j ) - (Kt 
*Deltat*ax)/DeltaSn) ; 

sigmaht->SetCellContent(i , j+1, sigmaht— >GetCellContent ( i , j +1) + 
(Kt*Deltat*ax)/DeltaSn) ; 

} 

} 

//*/ 

//Discharge currents: 
for (int j = 1; j <=ny ; j ++) { 
for (int i = 1; i <=nx ; i ++) { 

if (connh->GetCellContent(i , j ) < le20) { 

Rcon = ( 1/(2 * sigmas ) ) + connh— >GetCellContent(i , j ) ; 
It = ((-l)*(getV(l + nx*(j -1)) - 0)) / Rcon; 

sigmaht— >SetCellContent(i , j, sigmaht— >GetCellContent ( i , j ) + (I 
t*Deltat) /DeltaSn) ; 

//std::cout « "It = " « It « std : : endl ; 

} 

} 

} 

//Update sigmavt: 

for (int isigma = 1; isigma <= Nm; isigma++){ 

TString fname = TString :: Format ( "%sVsigma%i " , dir , isigma); 

TFile f (fname, "UPDATE"); // open the output file 

char labelsigma [ 1 0] ; 

sprintf ( labelsigma , " sigma%i " , isigma); 

TVectorD * sigmav ; 

f . GetObject( labelsigma , sigmav); 

for (int j = 1 ; j <= ny ; j ++) { 
for (int i = 1; i <= nx ; i++){ 

dummy = ((j — l)*nx+i — l)/nsub + 1; 

if (dummy == isigma) { 

(* sigmav ) ((j — l)*nx+i — ( isigma — 1)* nsub ) = sigmaht— >GetCellCo 
ntent(i, j); 

} 
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} 

} 



sigmav— >Write ( labelsigma ) ; 
f .CloseO ; 

} 

} 

TFile f3(fnameh, "UPDATE"); // open the output file for histograms 
char sigmahfName [ 1 6] ; 

s p r i n t f ( sigmahfName , " sigmah%i " , t S t ep i + nof S tep s ) ; 
TH2 *sigmahf; // Final surface charge density histogram 

sigmahf = new TH2D( Form (" Sigmah%i (x , y) t S tep i + nof Steps ) , "Surface 
charge density", nx , — lx/2, lx/2, ny , — ly/2, ly/2) ; 

//Fill the final surface charge density histogram: 
for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 
sigmahf — >SetCellContent(i , j, sigmaht— >GetCellContent ( i , j ) ) ; 

} 

} 

sigmahf —>Write ( sigmahfName) ; // Write the histogram after transient 
calculation 

delete sigmahf; 

f3 . Close ; 

delete sigmaht; 
delete connh; 

} 

#include " . . / Chani .h" 
#include "../ classes /Readp . h" 

void transient2 ( int tStepi , int nofSteps , double Deltat){ 
Readp rp ; 

double Ix = rp.GetlxO; 

double ly = rp.GetlyO; 

double ax = rp.Getax(); 

double ay = rp.GetayO; 

int nx = rp . Getnx ; 

int ny = rp.GetnyO; 

int Nm = rp.GetNm() ; 

double sigmas = rp . Getsigmas () ; 

int N = nx*ny; 

double DeltaSn = ax*ay; 

//stdi.cout « "Reading of the I— matrix" « std::endl; 
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// Read the I— matrix : 
TMatrixD 1 (1 ,N,1 ,N) ; 
int nsub = N / Nm; 

for (int i = 1 ; i <= Nm; i ++) { 
for (int j = 1; j <= Nm; j ++) { 

TString fnamel = TString :: Format ( "%slmatrix%i . root ", dir , i+(j-l)* 
Mn) ; 

TFile f (fnamel); // open the input file. 
char label [10]; 

sprintf ( label , " l-matrix%i " , i+(j-l)*Nm); 

TMatrixD *lp; // pointer of the sub matrix 
f . GetObject(label , Ip ) ; 

for (int m = ( j — l)*nsub + 1 ; m <= j*nsub; m++) { 
for (int n = (i — l)*nsub + l; n<= i*nsub; n++){ 

int mm = m— (j— l)*nsub; 
int nm = n — (i— l)*nsub; 

l(m, n) = (* Ip ) (mm,nm) ; 

// std : : cout « « m « " ," « n « ") = " « l(m,n) « " " 



] 

] 

delete Ip ; 

f .CloseO ; 

} 

} 

// std :: cout « "Reading of the I— matrix is completed" « std : : endl ; 



//Read the initial surface charge density histogram: 
TH2 *sigmaht; // Surface charge density histogram 

sigmaht = new TH2D( " SigmahtStep (x , y ) " , "Surface charge density" , n 
X, -lx/2, lx/2, ny, -ly/2, ly/2) ; 

TH2 *connh; //The connector histogram : 

connh = new TH2D( " Connectors " , "Connector histogram" , nx , —lx/2, 1 
x/2, ny, -ly/2, ly/2) ; 

cliar fnameh[100]; 

sprintf (fnameh , "%sh.root", dir); 

TFile f (fnameh, "UPDATE"); // open the output file for histograms 

TH2 *sigmahtd ; 
TH2 *connhd; 
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char objName [ 1 6 | ; 

sprintf (objName , "sigmah%i", tStepi); 
///. GetListOfKeys ()->Print() ; 

f . GetObject (objName , sigmahtd); // Read the initial surface charge 
density histogram 
f . GetObject (" connh" , connhd) ; // Read the connector histogram 



TVectorD sigmavt ( 1 ,N) ; // sigma vector 
TVectorD Vvt(l,N); // potential vector 

//Fill the sigma vector: 



for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 
sigmaht— >SetCellContent(i , j , 

} 

} 



sigmahtd ->GetCellContent ( i , j ) ) ; 



for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 
connh— >SetCellContent(i , j, connhd— >GetCellContent(i , j ) ) 

} 

} 

delete sigmahtd ; 
delete connhd; 

for (int j = 1; j <=ny ; j ++) { 
for (int i = 1; i<=nx; i++){ 
sigmavt( i+nx*(j — 1)) = sigmaht— >GetCellContent ( i , j ) ; 

} 

} 

double Kt ; // Surface current density variable 
double It ; // Discharge current variable 
double Rcon; // Connector resistance variable 
double sumt; // Charge sum variable 



for (int it = tStepi + 1; it <= tStepi + nofSteps ; it++){ 



std::cout « "time step " « it « std : : end! ; 



Vvt = 1 * sigmavt; 



//Update the surface charge density histogram due to the currents : 



// Handling the currents in x: 
for (int j = 1; j <=ny ; j ++) { 
for (int i = 1; i<=nx — 1; i++){ 



Kt = (( - 1) *( Vvt( i+l + nx*(j -1)) - Vvt ( i+nx *(j — 1)) ) * sigmas ) / ax ; 
sigmaht->SetCellContent(i , j, sigmaht->GetCellContent ( i , j ) - (Kt 
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* Deltat *ay ) / DeltaSn ) ; 

sigmaht— >SetCellContent ( i +1 , j, sigmaht— >GetCellContent ( i +1 , j ) + 
(Kt*Deltat*ay)/DeltaSn) ; 

} 

} 

// Handling the currents in y: 
for (int j = 1; j<=ny-l; j++){ 
for (int i = 1; i<=nx; i++){ 

Kt = (( — 1) *( Vvt( i+nx*j ) — Vvt ( i+nx *(j — 1)) ) * sigmas ) / ay ; 
sigmaht->SetCellContent(i , j, sigmaht->GetCellContent ( i , j ) - (Kt 
*Deltat*ax)/DeltaSn) ; 

sigmaht->SetCellContent(i , j+1, sigmaht->GetCellContent ( i , j + 1) + 
(Kt*Deltat*ax)/DeltaSn) ; 

} 

} 

//Discharge currents : 
for (int j = 1; j<=ny; j++){ 
for (int i = 1; i <=nx ; i++){ 

if (connh->GetCellContent(i , j ) < le20) { 

Rcon = ( 1/(2* sigmas ) ) + connh— >GetCellContent(i , j ) ; 

It = ((-l)*(Vvt(l + nx*(j -1)) - 0)) / Rcon; 

//cout « "Rcon = " « Rcon « endl ; 
//cout « "It = " « It « endl; 

sigmaht->SetCellContent(i , j, sigmaht->GetCellContent ( i , j ) + (I 
t*Deltat) /DeltaSn) ; 

} 

} 

} 

//Update the sigmavt 
for (int j = 1; j<=ny; j++){ 
for (int i = 1; i <=nx ; i++){ 
sigmavt( i+nx*(j — 1)) = sigmaht— >GetCellContent(i , j ) ; 

} 

} 

// sumt = sigmavt . Sum( ) *DeltaSn ; 

//cout « "sum(" « it « ") = " « sumt « endl; // Sanity check 
(Charge Conservation). 

} 

char sigmahf Name [16]; 

s p r i n t f ( sigmahfName , "sigmah%i", t S t ep i + nof S tep s ) ; 
TH2 *sigmahf; // Final surface charge density histogram 

sigmahf = new TH2D(Form( " Sigmah%i (x , y) ", tS tepi + nof Steps ) , "Surface 
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charge density", nx , -lx/2, lx/2, ny , -ly/2, ly/2) ; 

//Fill the final surface charge density histogram : 
for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 

sigmahf — >SetCellContent(i , j, sigmaht— >GetCellContent ( i , j ) ) ; 

//std::cout « sigmahtd—>GetCellContent{i,j) « std : : endl ; 

} 

] 

sigmahf —>Write ( sigmahfName) ; // Write the histogram after transient 
calculation 

delete sigmahf; 

f. Close ; 

delete sigmaht; 
delete connh; 

} 

#ifndef _GETVC_ 
# define __GETVC_ 
//#include "../Chani.h" 
#include "../ classes /Readp . h" 

// This function reads and returns the potential of the i—th cell 
double getV(int i){ 
Readp rp ; 

int nx = rp . Getnx ; 
int ny = rp.GetnyO; 
int Nm = rp.GetNm(); 
double sigmas = rp . Getsigmas () ; 

int N = nx*ny; 
int nsub = N/Nm; 

int iV = ( i — l)/nsub + 1; //indice of the v sub— vector 
double V; // Potential value to return 

TString fname = TString :: Format ( "%sVsigma%i " , dir , iV) ; 

TFile f (fname); // open the output file of vectors 

TVectorD *Vp; // pointer for the V— vector 

char labelV[5]; 

sprintf (labelV , "Wei", iV) ; 

f . GetObject(labelV , Vp) ; 

V = (*Vp)(i - (iV-1) * nsub); 
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delete Vp; 

f .CloseO ; 

return V; 

} 

#endif 

#include " . . / Chani .h" 
#include "../ classes /Readp .h" 

double getQtotal ( int tStep){ 

Readp rp ; 

double ax = rp.Getax(); 
double ay = rp.GetayO; 
double nx = rp.Getnx(); 
double ny = rp . Getny () ; 

double sum = 0; 

TString fnameh = TString :: Format ( "%sh . root " , dir); 
TFile f (fnameh); // open the output file for histograms 
TH2 *sigmah; // Surface charge density histogram 
char objName[161; 

sprintf (objName , " sigmah%i " , tStep); 

f . GetObject (objName , sigmah); // read the surface charge distributi 
on at time = Deltat * tStep 

for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 
sum = sum + sigmah — >GetCellContent ( i , j ) ; 

} 

} 

delete sigmah ; 

sum = sum * ax * ay; //from charge density to total charge 
f. Close ; 
return sum; 

} 

#include " . ./Chani .h" 
#include " . . / globalConsts . h" 
#include "../ classes /Readp . h" 

double capacitance () { 

double C = 0; 
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Readp rp ; 

double Ix = rp.GetlxO; 
double ly = rp.GetlyO; 
double ax = rp.Getax(); 
double ay = rp.GetayO; 
int nx = rp.GetnxO; 
int ny = rp.GetnyO; 
int Nm = rp .GetNm() ; 

int N = nx*ny ; 

double DeltaSn = ax*ay; 

std : : cout « "Reading the 1— matrix" « std : : endl ; 

TMatrixD 1(1 ,N,1 ,N) ; 
int nsub = N / Nm; 

for (int i = 1 ; i <= Nm; i++){ 
for (int j = 1; j <= Nm; j++){ 

TString fnamel = TString :: Format ( "%slmatrix%i . root ", dir , i+(j— 1)* 
Nm) ; 

TFile f (fnamel); // open the input file. 
char label[12]; 

sprintf ( label , " l-matrix%i " , i+(j-l)*Nm); 

TMatrixD *lp; // pointer of the sub matrix 
f . GetObject(label , Ip) ; 

for (int m = (j — 1)* nsub + 1 ; m <= j*nsub; m++) { 
for (int n = (i — l)*nsub + l; n<= i*nsub; n++){ 

int mm = m— (j— l)*nsub; 
int nm = n — (i— l)*nsub; 

l(m, n) = (* Ip ) (mm,nm) ; 

} 

} 

delete Ip ; 

f .CloseO ; 

} 

} 

std :: cout « "Reading of the 1— matrix is completed" « std :: endl ; 
TMatrixD linv (1 ,N,1 ,N) ; 
linv = 1 ; 

std :: cout « "Calculating the inverse matrix" « std :: endl ; 
linv . Invert () ; 
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std::cout « "Matrix inversions is completed" « std : : endl ; 

TVectorD V(l ,N) ; 
TVectorD sigma(l,N); 

V += 1 ; // su rfa ce with unit potential 

sigma = linv * V; // surface charge density distribution 
C = sigma. Sum() * DeltaSn ; // capacitance 
TH2D * sigmah ; 

sigmah = new TH2D( " sigma (x , y) " , "Charge distribution corresponding 
to unit potential", nx , -lx/2, lx/2, ny , -ly/2, ly/2); 

for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i <=nx ; i++){ 
sigmah— >SetCellContent(i , j, sigma ( i+nx *(j —1)) ) ; 

} 

} 

TString fnameh = TString :: Format ( "%sh . root " , dir); 

TFile f (fnameh, "UPDATE"); // open the output file for histograms 

sigmah —>Write (" sigmahVconst ") ; // Write the histogram after transie 
nt calculation 

f. Close ; 

delete sigmah; 

return C; 

} 

#include " . . / Chani .h" 
#include " . . / globalConsts . h" 
#include "../ classes / Findij . h" 

// This function returns the total charge on a space point x,y at th 

e 

// time step tStep 

double getQxy ( double x, double y, int tStep){ 

double Qxy; 

Findij fijQxy(x,y); 
int i = fij Qxy . Geti () ; 
int j = fijQxy . Getj () ; 
double ax = fijQxy . Getax () ; 
double ay = fijQxy . Getay () ; 

TString fnameh = TString :: Format ( "%sh . root " , dir); 
TFile f (fnameh); // open the output file for histograms 
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TH2 Ksigmah; // Surface charge density histogram 
char objName[16]; 

sprintf (objName , " sigmah%i " , tStep); 

f . GetObject(objName , sigmah); // read the surface charge distributio 
n at time = Deltat * tStep 

Qxy = sigmah— >GetCellContent ( i , j ) ; 

delete sigmah; 

f .CloseO ; 

return Qxy; 

} 

#ifndef _READPH_ 
# define _READPH_ 

#include " . . / Chani .h" 

class Readp{ 

public : 

double Ix , ly , ax , ay , sigmas ; 
int nx , ny , Nhi; 
Readp ( ) ; 

double Getlx () ; 
double Getly () ; 
double Getax () ; 
double Getay () ; 
int Getnx () ; 
int Getny () ; 
int GetNmO ; 
double GetsigmasO; 

}; 

#endif // __READPH__ 

#ifndef _READPC_ 
# define _READPC_ 

#include "../ Chani. h" 
#include "Readp. h" 

Readp : : Readp () { 
char fname[100]; 

sprintf (fname , "%sparameters . root " , dir); 

TFile f (fname); // open the paramsxy . root file. EDIT THIS LINE ACC 
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ORDING TO YOUR WORKING DIRECrORY . 

// f . GetListOfKeys ( )—>Print ( ) ; // Print to see if the objects are r 
etrieved properly 

THl *params ; // The histogram to hold Ix , nx , ax, ly , ny , ay value 
s in 

f. GetObject( "params " ,params ) ; 
Ix = params —>GetBinContent ( 1 ) ; 
nx = params —>GetBinContent (2) ; 
ax = params— >GetBinContent (3) ; 
ly = params —>GetBinContent (4) ; 
ny = params— >GetBinContent (5) ; 
ay = params —>GetBinContent (6) ; 
Ntn = params —>GetBinContent (7) ; 
sigmas = params— >GetBinContent (8) ; 

f. Close ; 

} 

double Readp : : Getlx { return Ix ; } 

double Readp :: Getly { return ly ; } 

double Readp :: Getax { return ax; } 

double Readp :: Getay { return ay; } 

int Readp :: Getnx { return nx ; } 

int Readp :: Getny { return ny ; } 

int Readp :: GetNm() { return Ntn; } 

double Readp :: Getsigmas { return sigmas; } 

#endif // _READPC__ 

#ifndef __F1NDIJH__ 
#define __FINDIJH__ 

#include "Readp. h" 
#include " . ./ Chani .h" 

class Findij : public Readp{ 
public : 
double X, y; 
int imin , jmin ; 

Findij ( ); 

Findij ( double , double); 

int Geti () ; 
int Getj () ; 

}; 

#endif // __FINDIJH__ 

#ifndef __F1ND1JC__ 
#define __F1ND1JC__ 

#include " ../ Chani. h" 
#include "Findij .h" 



75 



Findij :: Findij ( double x, double y){ 
//Reading x & y vectors: 
char fnamev [ 1 00] ; 

sprintf (fnamev , "%sxy.root", dir); 

TFile f (fnamev); // open the paramsxy . root file. EDIT THIS LINE ACC 
ORDING TO YOUR WORKING DIRECTORY . 

TVectorD *xp; // Pointer for x— vector 
f . GetObject("x-vector" ,xp) ; 

TVectorD *yp; // Pointer for y~vector 
f . GetObject("y-vector" ,yp) ; 

int N = nx*ny ; 

TVectorD xv(l ,N) ; 
XV = *xp; 
delete xp ; 
TVectorD yv(l ,N) ; 
yv = *yp; 
delete yp; 

f .CloseO ; 

/ / Dummy variables : 

double dist , distmin ; 

for (int j = 1; j < ny+1; j ++) { 
for (int i = l; i < nx + 1; i++){ 

dist = pow ( ( XV ( i +nx * ( j — 1) )— x) ,2) + pow ( ( yv ( i +nx * ( j — 1) )— y ) ,2) ; 

//cout « "dist = " « dist « " i, j = " « i « ", " « j « " x 

V, yv = " « xv( i+nx*(j -1)) « " , " « yv( i+nx*(j -1)) « endl ; 

if (i == 1 && j == l){distmin = dist; imin = i; jmin = j;} 
if (dist <= distmin ){ distmin = dist; imin = i; jmin = j;} 

// please note that , 

// if the connector coordinates are in the middle of two cells , th 
e one comes late is automatically choosen 

} 

} 

//cout « "imin, jmin = " « imin « ", " « jmin « endl; 
} 

int Findij :: Geti { return imin; } 
int Findij :: Getj { return jmin; } 

# end if // __FINDIJC__ 
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APPENDIX B Definition and simulation files of Sections 3.3 & 3.4 



#include <iostream> 
#include "Chani.h" 

/* 

Dimensions and number of divisions are defined in this file 
All lengths are in meters 

Once the macro is run, I— matrix is calculated 
//*/ 

void definitions { 

double Ix = 2e— 2; // length in x direction 
int nx = 15; // number of divisions in x 
double ly = 2e— 2; // length in y direction 
int ny = 15; // number of divisions in y 

double sigmas = le— 5; // surface conductivity 

int Nm = 1 ; // I— matrix is calculated in Nm x Nm part 

parameters ( Ix , nx , ly , ny , Nm, sigmas); // call the function for par 
ameter definitions 

computelmn ; // call the function that computes the I— matrix 
} 

#include <iostream> 

#include "Chani.h" 

#include "THl" 

#include "TH2" 

#include " classes /Readp . h" 

void connectors(){ 



Readp rp ; 



double 


Ix 




rp 


Getlx 


double 


ly 




rp 


. Getly 


double 


ax 




rp 


. Getax 


double 


ay 




rp 


Getay () 


double 


nx 




rp 


. Getnx 


double 


ny 




rp 


. Getny 


double 


y = 




-ly/2 + ay/2; 


while 


(y 


< 


le 


-2){ 



addConnector(— Ix /2 , y, 0) ; 
y = y + ay; 
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} 



} 

#include <iostream> 
#include "Chani.h" 
#include "THl" 
#include "TH2" 

void initialcharge () { 

Readp rp ; 

double Ix = rp.GetlxO; 
double ly = rp.GetlyO; 
double ax = rp.Getax(); 
double ay = rp.GetayO; 
int nx = rp.GetnxO; 
int ny = rp.Getny(); 
int Mm = rp.GetNm(); 

int N = nx*ny; 

double DeltaSn = ax*ay; 

double qadd = ( 1 e4 * 1 .6e — 19) ; //add le4 elementary charge 
//5 areas due to the 2d gaussian probabilities : 

double aO = (1 - ROOT: : Math :: normal_cdf (-ax/2 , ax, 0)*2)*(1 - ROOT:: 
Math : : normal_cdf(-ay /2 , ay, 0)*2); 

double al = ROOT: : Math :: normal_cdf_c(ax/2 , ax ,0) *(1 - ROOT: : Math :: nor 

mal_cdf( — ay /2 , ay, 0)*2); 

double a2 = ROOT: : Math :: normal_cdf_c ( ax /2 , ax , 0) *ROOT: : Math :: normal_c 
df_c(ay/2,ay ,0) ; 

double a3 = ROOT: : Math :: normal_cdf_c(ay/2 , ay ,0) *(1 - ROOT: : Math :: nor 
mal_cdf(— ax/2 , ax, 0)*2); 

//Add the initial charge with the gaussian probabilities: 

addCharge(0, 0, qadd*aO, 0); // C 

addCharge(ax , 0, qadd*al , 0); // R 

addCharge ( ax , ay, qadd*a2, 0); // UR 

addCharge(0, ay, qadd*a3 , 0); // U 

addCharge(— ax , ay, qadd*a2 , 0); // UL 

addCharge(— ax , 0, qadd*al , 0); // L 

addCharge(— ax , —ay, qadd*a2 , 0); // DL 

addCharge(0, -ay, qadd*a3 , 0); // D 

addCharge(ax , -ay, qadd*a2 , 0); // UR 

] 
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APPENDIX B.l Simulation files of Section 3.3 



#include <iostream> 
#include "Chani.h" 
#include "THl" 
#include "TH2" 
#include " connectors .C" 



void main { 



for (int i = 0; i < 2000; i++){ 
transient2 (i , 1 , le-10) ; 

} 



#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 



<iostream > 

"Chani .h" 

"THl" 

"TH2" 

" TCanvas " 

" TMath" 

"TMatrixD" 

"TVectorD" 

"math.h" 

"TFile" 



void getresults () { 



Readp rp ; 

double Ix = rp.GetlxO 

double ly = rp.GetlyO 

double ax = rp.Getax() 

double ay = rp.GetayO 

int nx = rp.GetnxO; 

int ny = rp.GetnyO; 



cout « " ax = " « ax « endl ; 



double Deltat = le-11; 



int tStepl = 0; 

int tl = tStepI * Deltat * le9; 

int tStep2 = 200; 

int t2 = tStep2 * Deltat * le9; 

int tStep3 = 500; 

int t3 = tStep3 * Deltat * le9; 

int tStep4 = 800; 

int t4 = tStep4 * Deltat * le9; 



TH2 *sigmahdl; // Surface charge density histogram for drawing 
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sigmahdl = new TH2D( "SigmahO(x , y) " , Form( " Surf ace charge density at 
t = %i ns", tl) , nx, -lx/2, lx/2, ny , -ly/2, ly/2) ; 

TH2 *sigmahd2; // Surface charge density histogram for drawing 
sigmahd2 = new TH2D( "SigmahlO(x , y) " , Form( " Surface charge density at 
t = %i ns", t2) , nx, -lx/2, lx/2, ny , -ly/2, ly/2) ; 

TH2 *sigmahd3; // Surface charge density histogram for drawing 
sigmahdS = new TH2D( "SigmahlOO(x , y) " , Form( " Surf ace charge density a 
t t = %i ns", t3) , nx, -lx/2, lx/2, ny , -ly/2, ly/2) ; 

TH2 *sigmahd4; // Surface charge density histogram for drawing 
sigmahd4 = new TH2D( "Sigmah200(x , y) " , Form( " Surf ace charge density a 
t t = %i ns", t4) , nx, -lx/2, lx/2, ny , -ly/2, ly/2) ; 

char fnameh [ 1 00] ; 

sprintf (fnameh , "%sh.root", dir); 

TFile f (fnameh); // open the output file for histograms 

f . GetListOfKeys ()— >Print ; // Print to see if the objects are retri 

eved properly 



TH2 


* sigmahl ; 


// 


S u rfa c e 


charge 


density 


histogram 


TH2 


* sigmah2 ; 


// 


Surface 


charge 


density 


histogram 


TH2 


* sigmahS ; 


// 


Surface 


charge 


density 


histogram 


TH2 


*sigmah4 ; 


// 


Surface 


charge 


density 


histogram 



char objNamel[16]; 

sprintf (objNamel , " sigmah%i " , tStepl); 
char objName2 [ 1 6 ] ; 

sprintf (objName2 , " sigmah%i " , tStep2); 
char objName3 [ 1 6 ] ; 

sprintf (objNameB , " sigmah%i " , tStep3 ) ; 
char objName4 [ 1 6 ] ; 

sprintf (objName4 , " sigmah%i " , tStep4); 

f . GetObject(objNamel , sigmahl); // read the surface charge distribut 
ion at time = Deltat * tStep 

f . GetObject(objName2 , sigmah2); // read the surface charge distribut 
ion at time = Deltat * tStep 

f . GetObject(objName3 , sigmah3 ) ; // read the surface charge distribut 
ion at time = Deltat * tStep 

f . GetObject (objName4 , sigmah4); // read the surface charge distribut 
ion at time = Deltat * tStep 

//Fill the surface charge density histograms: 
for (int j = 1; j<=ny; j++){ 
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sigmahl ->GetCellContent(i , j ) ) ; 



for (int i = 1; i<=nx; i++){ 

sigmahdl— >SetCellContent(i , j 
} 

} 

for (int j = 1; j <=ny ; j ++) { 
for (int i = 1; i<=nx; i++){ 
sigmahd2->SetCellContent(i , j, sigmah2 ->GetCellContent ( i , j ) ) ; 

} 

} 

for (int j = 1; j <=ny ; j ++) { 
for (int i = 1; i<=nx; i++){ 
sigmahdS— >SetCellContent(i , j, sigmah3 — >GetCellContent(i , j ) ) ; 

} 



for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 
sigmahd4->SetCellContent(i , j 

} 

} 



sigmah4->GetCellContent(i , j ) ) ; 



f. Close ; 



sigmahdl->SetStats (0) 
sigmahd2->SetStats (0) 
sigmahd3->SetStats (0) 
sigmahd4->SetStats (0) 



TCanvas *cl = new TCanvas( "cl ", "cl " ,1000,600) ; 

cl->Divide (2 ,2) ; 
cl->cd(l) ; 

sigmahdl— >Draw( " colz " ) ; 



cl->cd(2) ; 

sigmahd2->Draw( " colz " ) ; 
cl->cd(3) ; 

sigmahdB— >Draw( " colz " ) ; 
cl->cd(4) ; 

sigmahd4->Draw( " colz " ) ; 



cl->cd() ; 



cout « t4 « endl ; 



} 

#include <iostream> 

#include "Chani.h" 

#include "THl" 

#include "TH2" 

#include "TCanvas" 

#include "TMath" 

#include "TMatrixD" 

#include "TVectorD" 

#include "math.h" 

#include "TFile" 
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void totalcharge { 
Readp rp ; 

double Ix = rp.GetlxO; 
double ly = rp.GetlyO; 
double ax = rp.Getax(); 
double ay = rp . Getay () ; 
int nx = rp.GetnxO; 
int ny = rp . Getny () ; 
int Mm = rp.GetNm(); 

//Total charge histogram: 

int tStepi = 0; // first time step 

int tStepf = 2000; // final time step 

double Deltat = le — 10; // duration of time steps 

int DeltatStep = 1 ; // histograms corresponding to 0, DeltatStep , 2* 
DeltatStep , 3*DeltatStep ... instances are recorded 

int nBinTotal = ( tStepf — tStepi )/ DeltatStep + 1; //total number of bi 
ns 

double ti = tS tepi * Deltat ; 
double tf = tStepf *Deltat ; 

THl * totalcharge ; 

totalcharge =newTHlD("Qvs t"," Total charge versus time ", nBinTotal 
, ti , tf ) ; 

char fnameh [ 1 00] ; 

sprintf (fnameh , "%sh.root", dir); 

TFile f (fnameh); // open the output file for histograms 

f . GetListOfKeys — >Print ; // Print to see if the objects are retri 
eved properly 

for (int it = 0; it <= tStepf; it += DeltatStep ){ 

TH2 *sigmah; // Surface charge density histogram 
char objName[16]; 

sprintf (objName , "sigmah%i", it); 

f . GetObject(objName , sigmah); // read the surface charge distributi 
on at time = Deltat * tStep 

double sum = 0; 

for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 
sum = sum + sigmah— >GetCellContent(i , j ) ; 

} 

} 
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sum = sum * ax * ay; //from surface charge density to total charge 
totalcharge — >SetBinContent (it/DeltatStep+1, sum) ; 
delete sigmah ; 

} 

f . CloseO ; 

totalcharge ->SetStats (0) ; 

totalcharge->SetXTitle("Time (s)") ; 

totalcharge->SetYTitle("Charge (Coulomb)") ; 

totalcharge —>Draw() ; 

//*/ 

} 
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APPENDIX B.2 Simulation file of Section 3.4.1 



#include <iostream> 

#include "Chani.h" 

#include "THl" 

#include "TH2" 

void CvsN(){ 

// Dimentions and initial number of divisons in x and y axes: 

double Ix = 2e— 2; 
double nxi = 7; 
double ly = 2e-2; 
double nyi = 7; 

//Relative Tolerance for Capacitance 
double tol = 0.001; 

//Absolut maximum number for increasing of nx & ny : 
int imax = 30; 

// nx/ny ratio should be held constant 

double ratio = nxi/nyi; // ly < Ix should be Choosen 

THl * Chd; 

Chd = new TH1D( "C EOVIMY" , "DLMMY Capacitance with increasing number 
of divisions", imax + 1, 0, imax); 

THl * nxhd; 

nxhd = new TH1I( "nx DOVIMY" , 'TOVIMY Increase in nx" , imax + 1, 0, imax) 



THl * nyhd; 

nyhd = new TH1I( "ny IXMVIY" , 'TOVIVIY Increase in ny" , imax + 1, 0, imax) 
THl * Nhd; 

Nhd = new TH11( "N DUMMY" , "DUMMY Increase in total number of divisio 
ns", imax+1, 0, imax); 

THl * Ehd; 

Ehd = new TH1D(" Error DOVIMY', "IXMMY Error vs increasing number of d 
ivisions", imax, 1, imax); 

int nx , ny ; 

double sigmas = le— 5; // surface conductivity ( irrelevant for this c 
alculation . ) 

double Nm = 1 ; // 1—Part matrix calculation 
double C, error ; 

for (int i = 0; i <= imax+1; i++) { 
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cout « "i = " « i « endl ; 



ny = nyi + i ; 

nx = nxi + ratio * i ; 

parameters ( Ix , nx , ly , ny , Nm, sigmas); // call the function for pa 
rameter definitions 

computelmnO ; // call the function that computes the I— matrix 

C = capacitance () ; 

Chd->SetBinContent (i , C) ; 
nxhd— >SetBinContent ( i , nx) ; 
nyhd— >SetBinContent ( i , ny) ; 
Nhd— >SetBinContent ( i , nx*ny) ; 

if (i >= 1){ 

error = fabs ((C - Chd->GetBinContent (i -1)) ) /Chd->GetBinContent (i -1 

); 

Ehd— >SetBinContent ( i , error); 



if (error < tol) { 

cout « "Tolerance for relative error is achieved for nx = " « n 
X « " ny = " « ny « endl ; 

break ; 

} 

} 

} 

THl * Ch; 

Ch = new TH1D("C vs N" , "Capacitance vs Number of divisions", i+1, n 
xi*nyi , nx*ny) ; 

for (int j = 0; j <= i ; j ++){ 
Ch->SetBinContent ( j +1 , Chd->GetBinContent ( j ) ) ; 
cout « "Ch("« j « ")= " « Ch->GetBinContent (j ) « endl; 

} 

delete Chd; 
THl * nxh ; 

nxh = new THlD("nx", "Increase in nx" , i+1, 0, i); 

for (int j = 0; j <= i ; j++){ 
nxh->SetBinContent ( j +1 , nxhd->GetBinContent ( j ) ) ; 

} 

delete nxhd; 
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THl * nyh; 

nyh = new THlI("ny", "Increase in ny" , i+1, 0, i); 

for (int j = 0; j <= i ; j++){ 
nyh->SetBinContent ( j +1 , nyhd->GetBinContent ( j ) ) ; 

} 

delete nyhd; 

THl * Nh; 

Nh = new TH1I("N", "Increase in total number of divisions", i+1, 

i); 

for (int j = 0; j <= i ; j++){ 

Nh->SetBinContent ( j +1 , Nlid->GetBinContent ( j ) ) ; 

} 

delete Nhd; 
THl * Ehd; 

Eh = new TH1D( " Error " , "Relative Error vs Number of divisions", 
nxi*nyi , nx*ny) ; 

for (int j = 1; j <= i; j++){ 
Eh->SetBinContent (j , Ehd->GetBinContent ( j ) ) ; 

} 

delete Ehd; 
Ch->SetStats (0) ; 

Ch->SetYTitle("Capacitance (F)") ; 
Ch->SetXTitle( "Number of cells"); 
Ch->SetTitleOffset (1.4 , "Y" ) ; 
nxh->SetStats (0) ; 
nyh->SetStats (0) ; 
Nh->SetStats (0) ; 
Eh->SetStats (0) ; 

Eh->SetXTitle( "Number of cells"); 

Eh— >SetYTitle (" Relative error in capacitance"); 

Eh->SetTitleOffset (1.4 , "Y" ) ; 

TCanvas *cl = new TCanvas( " cl " , " cl " ,600 ,400) ; 

cl->cd(l) ; 

Ch->Draw() ; 

cl->cd() ; 

TCanvas *c2 = new TCanvas( " c2 " , " c2 " ,600 ,400) ; 
c2->cd(l) ; 

Eh->Draw() ; 
c2->cd ; 
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TH2 *sigmah; // Surface charge density histogram for drawing 

sigmah = new TH2D( " Sigma (x , y) " , "Surface charge distribution corresp 

onding to unit potential" , nx , — lx/2, lx/2, ny , — ly/2, ly/2) ; 

TH2 *sigmahd; // Surface charge density histogram for drawing 

char fnameh[100]; 

sprintf ( fnameh , "%sh.root", dir); 

TFile f (fnameh); // open the output file for histograms 

f . GetObject ( " sigmah Vconst " , sigmahd ) ; 

for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 
sigmah— >SetCellContent(i , j, sigmahd— >GetCellContent(i , j )) ; 

} 

} 

f. Close ; 

sigmah->SetStats (0) ; 

TCanvas *c3 = new TCanvas( " c3 " , " c3 " ,600 ,600) ; 
c3->cd(l) ; 

sigmah— >Draw( " colz " ) ; 

c3->cd() ; 

} 
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APPENDIX B.3 Simulation file of Section 3.4.2 



#include <iostream> 
#include "Chani.h" 
#include "THl" 
#include "TH2" 
#include " connectors .C" 
#include " initialcharge .C" 

void tvsN { 

// Dimentions and initial number of divisons in x and y axes: 

double Ix = 2e— 2; 
double nxi = 7; 
double ly = 2e-2; 
double nyi = 7; 

//Relative Tolerance for Discharge Time 
double tol = -0.001; 

//Absolut maximum number for increasing of nx & ny : 
int imax = 20; 

// nx/ny ratio should be held constant 

double ratio = nxi/nyi; // ly < Ix should be Choosen 

THl * tdhd; 

tdhd = new TH1D( " t DUMVIY" , "DOVIIVIY Discharge time with increasing num 
ber of divisions", imax + 1, 0, imax); 

THl * nxhd ; 

nxhd = new TH1I( "nx DLMVIY" , "DUVDVIY Increase in nx" , imax + 1, 0, imax) 



THl * nyhd; 

nyhd = new TH1I( " ny DOVIVIY" , "D0V1MY Increase in ny" , imax + 1, 0, imax) 
THl * Nhd; 

Nhd = new TH1I( "N DUMMY" , "DOVMY Increase in total number of divisio 
ns", imax+1, 0, imax); 

THl * Ehd; 

Ehd = new TH1D(" Error DUMMY", "DUMMY Error vs increasing number of d 
ivisions", imax, 1, imax); 

int nx , ny , N; 

double qadd = le4*1.6e — 19; 

double sigmas = le— 5; // surface conductivity 
double Mm = 1 ; // 1—Part matrix calculation 
double C, error, tdischarge ; 
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double Deltat = le-10; 

for (int i = 0; i <= imax+1; i += 2) { 

cout « "i = " « i « endl ; 

ny = nyi + i ; 

nx = nxi + ratio * i ; 

parameters ( Ix , nx , ly , ny , Nm, sigmas); // call the function for pa 
rameter definitions 

computelmnO ; // call the function that computes the I— matrix 

connectors () ; 

addCharge(0, 0, qadd , 0); 

cout « "qadd = " « qadd « endl; 

for (int it = 0; it < 2000; it++){ 

transient2 ( it , 1, Deltat); 
double Q = getQtotal ( it +1) ; 

if (Q/qadd < 0.1) { break;} 

} 

tdischarge = Deltat*it; 

tdhd->SetBinContent (i/2+1, tdischarge); 
nxhd->SetBinContent ( i /2+1 , nx) ; 
nyhd->SetBinContent ( i /2+ 1 , ny) ; 
Nhd->SetBinContent(i/2+l , nx*ny) ; 

if (i/2 >= 1){ 

cout « "tdischarge = " « tdischarge « endl; 

cout « "tdhd( " « i/2 « ") = " « tdhd->GetBinContent ( i /2) « e 
ndl; 

error = (tdischarge - tdhd->GetBinContent ( i /2) ) / tdhd->GetBinCont 
ent(i/2) ; 

error = fabs(error); 

cout « " error = " « error « endl ; 
Ehd->SetBinContent ( i /2 , error); 

if (error < tol) { 
cout « "Tolerance for relative error is achieved for nx = " « n 
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X « " ny = " « ny « endl ; 
break ; 
} 

} 



} 

THl * tdh; 

tdh = new TH1D( " t_discharge vs N" , "Discharge time vs Number of divi 
sions", i/2+1, nxi*nyi , nx*ny) ; 

for (int j = 0; j <= ill; j++){ 
tdh->SetBinContent(j+l, tdhd->GetBinContent ( j +1)) ; 
cout « "tdh("« j « ")= " « tdh->GetBinContent ( j +1) « endl; 

} 

delete tdhd ; 

THl * nxh ; 

nxh = new THlD("nx", "Increase in nx" , i+1, 0, i); 

for (int j = 0; j <= i/2; j ++) { 
nxh->SetBinContent (j +1 , nxhd->GetBinContent (j +1)) ; 

} 

delete nxhd; 
THl * nyh ; 

nyh = new THlI("ny", "Increase in ny" , i+1, 0, i); 

for (int j = 0; j <= i/2; j ++) { 
nyh->SetBinContent (j +1 , nyhd->GetBinContent (j +1)) ; 

} 

delete nyhd; 
THl * Nh; 

Nh = new TH1I("N", "Increase in total number of divisions", i+1, 0, 

i); 

for (int j = 0; j <= i/2; j++){ 

Nh->SetBinContent ( j +1 , Nhd->GetBinContent ( j +1) ) ; 

} 

delete Nhd; 
THl * Ehd; 

Eh = new TH1D( " Error " , "Relative Error vs Number of divisions", i, 
nxi*nyi , nx*ny) ; 

for (int j = 1; j <= i/2; j++){ 
Eh->SetBinContent (j , Ehd->GetBinContent ( j ) ) ; 

} 
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delete Ehd; 



tdh->SetStats (0) ; 

tdh->SetYTitle (" Discharge time (s)"); 
tdh->SetXTitle( "Number of cells"); 
tdh->SetTitleOffset (1 .4 , "Y" ) ; 
nxh->SetStats (0) ; 
nyh->SetStats (0) ; 
Nh->SetStats (0) ; 
Eh->SetStats (0) ; 

Eh->SetXTitle( "Number of cells"); 

Eh->SetYTitle("Relative error in discharge time"); 
Eh->SetTitleOffset (1.4 , "Y" ) ; 

TCanvas *cl = new TCanvas ( " cl " , " cl " ,600 ,400) ; 

cl->cd(l) ; 

tdh->Draw() ; 

cl->cd() ; 

TCanvas *c2 = new TCanvas( "c2" , "c2" ,600,400) ; 

c2->cd(l) ; 

Eh->Draw() ; 

c2->cd ; 

//*/ 

} 
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APPENDIX C Definition and simulation files of Section 3.5 



#include <iostream> 
#include "Chani.h" 

/* 

Dimensions and number of divisions are defined in this file 
All lengths are in meters 

Once the macro is run, I— matrix is calculated 
//*/ 

void definitions { 

double Ix = 4e— 2; // length in x direction 
int nx = 99; // number of divisions in x 
double ly = 4e— 2; // length in y direction 
int ny = 99; // number of divisions in y 

double sigmas = le— 5; // surface conductivity 

int Nm = 1 ; // I— matrix is calculated in Nm x Nm part 

parameters ( Ix , nx , ly , ny , Nm, sigmas); // call the function for par 
ameter definitions 

computelmn ; // call the function that computes the I— matrix 
} 

// this file adds a periodic charge distribution on the surface 

// sigma = sigmaO * sin(kx) e^(—t/tau) 

// Chani will simulate the relaxation of the charge. 

#include <iostream> 
#include "Chani.h" 
#include "TH2.h' 
#include "TMath.h" 

// compile and run with: 

// g++ 'root — config — glibs ' —I 'root — config — incdir ' initialcharge 
sin .C ; . / a . out 

//void initialchar ge sin ( ) { 
int initialchargesin { 

Readp rp ; 

double Ix = rp.GetlxO; 

double ly = rp.GetlyO; 

double ax = rp.Getax(); 

double ay = rp.GetayO; 

int nx = rp.GetnxO; 

int ny = rp.Getny(); 
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int Nm = rp.GetNm(); 

int N = nx*ny ; 

double DeltaSn = ax*ay; 

TString fnamev = TString :: Format ( "%sxy . root " , dir); 

// Reading of x & y vectors: 
TFile f (fnamev); 

TVectorD *xp; // Pointer for x— vector 
f . GetObject( "x-vector " ,xp) ; 
TVectorD x(l ,N) ; 
X = *xp; 
delete xp ; 

TVectorD *yp; // Pointer for y— vector 
f . GetObject( "y-vector " ,yp) ; 
TVectorD y(l ,N) ; 
y = *yp; 
delete yp; 

f .CloseO ; 

double k = 1000; 

cout « "k = " « k « endl ; 

double sigmaaddO = (1 e4 * 1 .6e — 19)/ DeltaSn ; //le4 elementary charge 
// add the initial periodic charge: 
double sigmaadd ; 
TH2 * sigmahO ; 

sigmahO = new TH2D( " sigmaO (x , y ) " , "Surface charge density at t = 0" , 
nx, -lx/2, lx/2, ny, -ly/2, ly/2) ; 

for (int j = 1; j <= ny; j++){ 
for (int i = 1 ; i <= nx ; i++){ 

//std::cout « "i = " « i « std::endl; 

sigmaadd = sigmaaddO*cos (k*x( i +(j — l)*nx) ) ; 
sigmahO— >SetCellContent(i , j , sigmaadd) ; 

} 

} 

TString fnameh = TString :: Format ( "%sh . root " , dir); 

TFile f 2 ( fnameh , "UPDATE" ) ; // open the output file for histograms 

sigmahO— >Write (" sigmahO ") ; // Write the surface with periodic initia 
1 charge 

f2 . Close ; 
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} 

#ifndef __CINT__ 
int main () { 

initialchargesin2 () ; 

return 0; 

} 

#endif 

#include <iostream> 
#include "Chani.h" 
#include "THl.h" 
#include "TH2.h" 
#include " connectors .C" 

int main { 

int tStepi = 0; 

int tStepf = 1000; 

int DeltatStep = 1; 

double Deltat = le-11; // 0.01 ns 

//double C = capacitance ; 

// std : : cout « "C = " « C « std : : endl ; 

for (int i = tStepi ; i < tStepf ; i = i + DeltatStep ){ 

transient2 (i , DeltatStep, Deltat); 

//cout « "Q total = " « getQtotal(i) « endl; 

} 

return 0; 

} 

#include <iostream> 
#include "Chani.h" 
#include "THl" 
#include "TH2" 
#include "TCanvas" 
#include "TMath" 
#include "TMatrixD" 
#include "TVectorD" 
#include "math.h" 
#include "TFile.h" 
#include " classes /Readp . h" 

void qxy(){ 

int tStepi = 0; // first time step 

int tStepf = 600; // final time step 

double Deltat = le — 11; // duration of time steps 

int DeltatStep = 10; // histograms corresponding to 0, DeltatStep , 
2*DeltatStep , 3*DeltatStep ... instances are recorded 

int nBinTotal = ( tStepf-tStepi )/ DeltatStep + 1; //total number of b 
ins 
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double ti = tStepi * Deltat ; 
double tf = tStepf *Deltat ; 



double X = 0; 
double y = 0; 

THl *qxy; 

qxy = new THIDC'Q vs t", "Charge at (0,0) versus time " , nBinTotal , ti 
*le9, tf*le9); 

double q; 

for (int it = tStepi ; it <= tStepf ; it += DeltatStep){ 
q = getQxy (x , y , it —tStepi ) ; 

//std::cout « "q at time step " « it « " is " « q « std : : endl 
qxy->SetBinContent ( it / DeltatStep + 1 , q) ; 

} 

qxy->SetStats (0) ; 
qxy->SetXTitle("Time (ns)"); 
qxy->SetYTitle( "Charge (Coulomb)") ; 

TCanvas *cl = new TCanvas( "cl " , "cl " ,600 ,400) ; 



qxy— >Draw () ; 
qxy->Fit ( " expo" ) ; 

TFl *fit = qxy->GetFunction("expo") ; 

//fit —>SetParameter ( , log ( qxy—>GetBinContent( tStepi ) ) ) ; 
// cl->Update() ; 
double tau ; 

tau = (-!)/( fit ->GetParameter(l) ) ; 
std::cout « "tau = " « tau « std :: endl ; 

} 
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APPENDIX D Definition and simulation files of Section 3.6 



#include <iostream> 
#include "Chani.h" 

/* 

Dimensions and number of divisions are defined in this file 
All lengths are in meters 

Once the macro is run, I— matrix is calculated 
//*/ 

void definitions { 

double Ix = lOe— 2; // length in x direction 
int nx = 667; // number of divisions in x 
double ly = 150e— 6; // length in y direction 
int ny = 1; // number of divisions in y 



double sigmas = 33.3e— 6; // surface conductivity 

int Mn = 1 ; // I— matrix is calculated in Nm x Nm part 

parameters ( Ix , nx , ly , ny , Nm, sigmas); // call the function for par 
ameter definitions 

computelmn ; // call the function that computes the I— matrix 
} 

#include <iostream> 
#include "Chani.h" 
#include "THl" 
#include "TH2" 

void connectors { 

addConnector(-5e-2, 0, 0); 

} 

#include <iostream> 
#include "Chani.h" 
#include "THl" 
#include "TH2" 

void initialcharge { 

Readp rp ; 

double Ix = rp.GetlxO; 
double ly = rp.GetlyO; 
double ax = rp.Getax(); 
double ay = rp.GetayO; 
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int nx = rp.GetnxO 
int ny = rp.GetnyO 
int Nm = rp.GetNm() 



int N = nx*ny ; 
double DeltaSn 



ax*ay ; 



double qadd = ( 1 e4 * 1 .6e — 19) ; //add le4 elementary charge 
//J areas due to the 2d gaussian probabilities : 
double aO = (1 - ROOT: : Math :: normal_cdf (-ax/2 , ax, 0)*2); 
double al = ROOT: : Math :: normal_cdf (-ax/2 , ax, 0); 

//Add the initial charge with the gaussian probabilities: 
addCharge(0, 0, qadd*aO, 0); // C 
addCharge(ax , 0, qadd*al , 0); // R 
addCharge(— ax , 0, qadd*al , 0); // L 

} 

#include <iostream> 
#include "Chani.h" 
#include "THl" 
#include "TH2" 



void main () { 
double C; 

for (int i = 0; i < 2000000; i += 1000){ 
transient2 (i , 1000, le-10); 

} 
} 



#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 



<iostream > 

"Chani .h" 

"THl" 

"TH2" 

"TCanvas" 

" TMath" 

"TMatrixD" 

"TVectorD " 

"math.h" 

"TFile" 



void getresults () { 
Readp rp ; 

double Ix = rp.GetlxO 
double ly = rp.GetlyO 
double ax = rp . Getax () 
double ay = rp.GetayO 
int nx = rp.GetnxO; 
int ny = rp.GetnyO; 
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cout « " ax = " « ax « endl ; 



double Deltat = le-10; 

int tStepl = 0; 

int tl = tStepI * Deltat * le9; 

int tStep2 = 10000; 

int t2 = tStep2 * Deltat * le9; 

int tStep3 = 1000000; 

int t3 = tStep3 * Deltat * le9; 

int tStep4 = 2000000; 

int t4 = tStep4 * Deltat * le9; 



TH2 *sigmahdl ; // Surface charge density histogram for drawing 
sigmahdl = new TH2D( " SigmahO(x , y ) " , Form( " Surface charge density at 
t = %i ns", tl) , nx, -lx/2, lx/2, ny , -ly/2, ly/2) ; 

TH2 *sigmahd2; // Surface charge density histogram for drawing 
sigmahd2 = new TH2D( " SigmahlO(x , y ) " , Form( " Surface charge density at 
t = %i ns", t2) , nx, -lx/2, lx/2, ny , -ly/2, ly/2) ; 

TH2 *sigmahd3; // Surface charge density histogram for drawing 
sigmahd3 = new TH2D( " SigmahlOO (x , y ) " , Form( " Surf ace charge density a 
t t = %i ns", t3) , nx, -lx/2, lx/2, ny , -ly/2, ly/2) ; 

TH2 *sigmahd4; // Surface charge density histogram for drawing 
sigmahd4 = new TH2D( " Sigmah200 (x , y ) " , Form( " Surface charge density a 
t t = %i ns", t4) , nx, -lx/2, lx/2, ny , -ly/2, ly/2) ; 

char fnameh[100]; 

sprintf (fnameh , "%sh.root", dir); 

TFile f (fnameh); // open the output file for histograms 

f . GetListOfKeys ()— >Print ; // Print to see if the objects are retri 

eved properly 



TH2 


* sigmahl ; 


// 


Surface 


charge 


density 


histogram 


TH2 


* sigmah2 ; 


// 


Surface 


charge 


density 


histogram 


TH2 


*sigmah3 ; 


// 


Surface 


charge 


density 


histogram 


TH2 


* sigmah4 ; 


// 


Surface 


charge 


density 


histogram 



cliar objNamel [ 1 6] ; 
sprintf (objNamel , 

cliar objName2 [ 1 6] ; 
sprintf (objName2 , 

cliar objName3 [ 1 6] ; 
sprintf (objName3 , 

char objName4 [ 1 6] ; 



sigmah%i", tStepI); 
sigmah%i", tStep2); 
sigmah%i", tStep3 ) ; 
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sprintf (objName4 , " sigmah%i " , tStep4); 



f . GetObject(objNamel , sigmahl); // read the surface charge distribut 
ion at time = Deltat * tStep 

f . GetObject(objName2 , sigmah2); // read the surface charge distribut 
ion at time = Deltat * tStep 

f . GetObject(objName3 , sigmah3 ) ; // read the surface charge distribut 
ion at time = Deltat * tStep 

f . GetObject(objName4 , sigmah4); // read the surface charge distribut 
ion at time = Deltat * tStep 

//Fill the surface charge density histograms : 
for (int j = 1; j<=ny; j++){ 
for (int 1=1; 1 <=nx ; l++){ 
sigmahdl— >SetCellContent(i , j, sigmahl — >GetCellContent ( i , j )) ; 

} 

} 

for (int j = 1; j<=ny; j++){ 
for (int i = 1; i<=nx; i++){ 
sigmahd2->SetCellContent(i , j, sigmah2->GetCellContent(i , j ) ) ; 

} 

} 

for (int j = 1; j<=ny; j ++) { 
for (int i = 1; i<=nx; i++){ 
sigmahd3— >SetCellContent ( i , j, sigmahS — >GetCellContent ( i , j ) ) ; 

} 

} 

for (int j = 1; j<=ny; j++){ 
for (int 1 = 1; 1 <=nx ; l++){ 
sigmahd4->SetCellContent( i , j, sigmah4->GetCellContent(i , j ) ) ; 

} 

} 



f .CloseO ; 



sigmahdl->SetStats (0) 
sigmahd2->SetStats (0) 
sigmahd3->SetStats (0) 
sigmahd4->SetStats (0) 



TCanvas *cl = new TCanvas (" cl "," cl ", 1000 ,600) ; 
cl->Divide (2,2) ; 
cl->cd(l) ; 

sigmahdl— >Draw( " colz " ) ; 



cl->cd(2) ; 

sigmahd2->Draw( " colz " ) ; 
cl->cd(3) ; 

sigmahd3— >Draw( " colz " ) ; 
cl->cd(4) ; 

sigmahd4->Draw( " colz " ) ; 
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cl->cd() 



cout « t4 « endl ; 



} 



#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 



<iostream > 

"Chani .h" 

"THl" 

"TH2" 

"TCanvas" 

" TMath" 

"TMatrixD" 

"TVectorD " 

"math.h" 

"TFile" 



void totalcharge () { 
Readp rp ; 

double Ix = rp.GetlxO; 
double ly = rp.GetlyO; 
double ax = rp.Getax(); 
double ay = rp.GetayO; 
int nx = rp.GetnxO; 
int ny = rp.GetnyO; 
int Nm = rp .GetNm() ; 

//Total charge histogram: 



int tStepi = 0; // first time step 

int tStepf = 400000; // final time step 

double Deltat = le — 10; // duration of time steps 

int DeltatStep = 1000; // histograms corresponding to 0, DeltatStep , 

2*DeltatStep , 3*DeltatStep ... instances are recorded 
int nBinTotal = ( tStepf — tStepi )/ DeltatStep + 1; //total number of bi 
ns 

double ti = tStepi * Deltat ; 
double tf = tStepf *Deltat ; 

THl * totalcharge ; 

totalcharge = new TH1D( "Q vs t", "Total charge versus time ", nBinTotal 
, ti , tf ) ; 

char fnameh [100]; 

sprintf (fnameh , "%sh.root", dir ) ; 

TFile f (fnameh); // open the output file for histograms 

f . GetListOfKeys ()— >Print ; // Print to see if the objects are retri 
eved properly 

for (int it = 0; it <= tStepf ; it += DeltatStep ){ 
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TH2 Ksigmah; // Surface charge density histogram 



char objName[16]; 



sprintf (objName , "sigmah%i", it); 

f . GetObject(objName , sigmah); // read the surface charge distributi 
on at time = Deltat * tStep 

double sum = 0; 

for (int j = 1; j <=ny ; j++){ 
for (int i = 1; i<=nx; i++){ 
sum = sum + sigmah— >GetCellContent(i , j ) ; 

} 

} 

sum = sum * ax * ay; //from surface charge density to total charge 
total charge — >SetBinContent (it/DeltatStep+1, sum) ; 
delete sigmah; 

} 

f. Close ; 

totalcharge->SetStats(0); 

totalcharge->SetXTitle("Time (s)") ; 

totalcharge ->S at YTitle (" Charge (Coulomb)") ; 

totalcharge — >Draw ; 

//*/ 

} 
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APPENDIX E Definition and simulation files of Section 3.7 

#include <iostream> 
#include "Chani.h" 



/* 

Dimensions and number of divisions are defined in this file 
All lengths are in meters 

Once the macro is run, I— matrix is calculated 
//*/ 



void definitions { 



double Ix = lOOe— 3; // length in x direction 

int nx = 279; // number of divisions in x 

double ly = 18e— 3; // length in y direction 

int ny = 51; // number of divisions in y 

double sigmas = le— 5; // surface conductivity 

int Ntn = 3; // I— matrix is calculated in Nm x Nm part 

parameters ( Ix , nx , ly , ny , Nm, sigmas); // call the function for par 
ameter definitions 

computelmnO ; // call the function that computes the I— matrix 
} 

#include <iostream> 
#include "Chani.h" 
#include "THl" 
#include "TH2" 

void connectors(){ 

double y = — 9e — 3; 

while (y < 10e-3){ 

addConnector(-500e-3, y, 0); 
y = y + 9e-3; 

} 



} 

#include <iostream> 
#include "Chani.h" 
#include "THl" 
#include "TH2" 



void initialcharge { 
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Readp rp ; 

double Ix = rp.GetlxO; 
double ly = rp.GetlyO; 
double ax = rp.Getax(); 
double ay = rp.GetayO; 
int nx = rp.GetnxO 
int ny = rp . Getny () 
int Nm = rp.GetNm() 

int N = nx*ny; 

double DeltaSn = ax*ay; 

double qadd = ( 1 e4 * 1 . 6 e — 19) ; //add 1 e4 elementary charge 
//5 areas due to the 2d gaussian probabilities: 

double aO = (1 - ROOT: : Math :: normal_cdf (-ax/2 , ax, 0)*2)*(1 - ROOT:: 

Math : : normal_cdf ( — ay /2 , ay, 0)*2); 

double al = ROOT: : Math :: normal_cdf_c(ax/2 , ax ,0) *(1 - ROOT: : Math :: nor 
mal_cdf(-ay/2 , ay, 0)*2); 

double a2 = ROOT: : Math :: normal_cdf_c(ax/2 , ax ,0) *ROOT: : Math :: normal_c 
df_c(ay/2,ay ,0) ; 

double a3 = ROOT: : Math :: normal_cdf_c(ay/2 , ay ,0) *(1 - ROOT: : Math :: nor 
mal_cdf(-ax/2 , ax, 0)*2); 

//Add the initial charge with the gaussian probabilities: 

addCharge(0, 0, qadd*aO, 0); // C 

addCharge(ax , 0, qadd*al , 0); // R 

addCharge(ax , ay, qadd*a2, 0); // UR 

addCharge(0, ay, qadd*a3 , 0); // U 

addCharge(— ax , ay , qadd*a2 , 0) ; // UL 

addCharge(— ax , 0, qadd*al , 0); // L 

addCharge(— ax , —ay, qadd*a2 , 0); // DL 

addCharge(0, -ay, qadd*a3 , 0); // D 

addCharge(ax , -ay, qadd*a2 , 0); // UR 

} 

#include <iostream> 
#include "Chani.h" 
#include "THl" 
#include "TH2" 

void main () { 

for (int i = 0; i < 200; i++){ 
transient(i, 1, le — 10); 

} 



} 
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